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ON THE COMPUTATION OF STRAIN AND 
DISPLACEMENT IN A PRISM PLASTICALLY STRAINED 
BY TORSION 
By R. V. SOUTHWELL (Ozford) 

[Received 20 September 1949] 


SUMMARY 

For an investigation (now in course of publication) of plastic distortion in some 
systems of plane strain and of plane stress, a physical hypothesis was needed to 
supplement the (shear stress)/(shear strain) diagram which is presumed in Prandtl’s 
treatment of plastic torsion. The chosen hypothesis was consistent with that 
diagram, being simplified in relation to torsion by the fact that in Prandtl’s solution 
the directions of principal shear stress do not alter once plastic straining has begun. 
Strains and displacements (not considered in his treatment) were found to be readily 
deducible from the torsional stresses. 

This conclusion is confirmed by the work of P. G. Hodge (1). But Hodge restricts 
his detailed treatment to ‘the fully plastic state’, whereas here the need thus to 
restrict is obviated by a use of Relaxation Methods for determination of the plastic- 
elastic interface. Detailed computations have been made by Miss G. Vaisey for a 
hollow rectangle treated by F.S. Shaw (10). The axial (warping) displacement is 
exhibited by contours, also the plastic strains within the enclaves (regions of plastic 
straining). 


Introduction 
1, L. PRANDTL (2) extended so as to take account of non-elastic straining 
Saint-Venant’s well-known theory (3) of torsion in an isotropic prism. 
His treatment presumes the relation between shear stress and strain to 
be of the kind which is depicted in Fig. 1. For shear stress less than some 
given value f;- the usual (Hooke’s) relation holds and in consequence Y’, 
the ‘stress-function’, satisfies the equation 
Vy+2=0 (1) 

as in Saint-Venant’s theory. But the shear stress cannot exceed this 
limiting value, with which, accordingly, any shear strain greater than yy 
can be associated; so an upper limit is imposed upon the gradient of ¥’. 

The work of Boussinesq and Filon (see (4), § 219) shows that the shear- 
stress found from (1) attains its greatest value on the cylindrical boundary. 
Therefore (in Prandtl’s argument) knowing Y on the boundary we can 
impose upper limits on its value at neighbouring points. They can be 
represented by a surface of constant slope, attached to the boundary in 
the way that a roof is attached to the wall of a house. (This surface is 
ommonly termed a ‘Prandtl roof’.) 


Earlier (1903) Prandtl (5) had propounded his ‘membrane analogue’ 


[Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 4 (1949)] 
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whereby ‘’, in (1), is visualized as the transverse displacement of a mem. 
brane stretched with uniform tension and fixed (so that ‘’ = 0) at the 
given boundary. This has been made the basis of experiments with soap- 
films to determine torsional stresses in elastic prisms (shafts) of shapes 
which are not tractable by mathematics;* also (in conjunction with his 
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Shear strain 
Fig. 1. 
notion of a ‘roof’) of experiments to determine the altered distribution 
which results from overstrain. 

2. Christopherson and Southwell (7) replaced (1) by an approximation 
in finite differences, and on that basis showed that torsional stresses 
whether plastic or elastic can be determined, for a shaft of any section, 
without recourse to experiment. Various shapes of cross-section have 
been thus treated, by Christopherson (8), Shaw (9, 10), Allen e¢ al. (11), 
and Cassie and Dobie (12) as regards elastic, also by Christopherson (8 
and by Shaw (10) as regards plastic straining. Shaw dealt with some 
hollow (multiply-connected) sections, and also discussed the computation 
of ‘residual’ stresses (those which persist after the applied torque has 
been removed). 

3. Thus methods for determining plastic stresses are now fairly widel) 
known. But Hodge (1) alone seems to have considered the accompanying 
strains (that is, the nature of the consequent distortion), and his detailed 
treatment is restricted to the ‘fully plastic state’ in which (§ 1) the ‘mem- 
brane’ is almost everywhere in contact with the ‘roof’ and only a thin core 
of material remains elastic.+ Here a slightly different treatment is pre- 
sented, devised (in ignorance of Hodge’s work) to supplement a longer 

* Cf. Southwell (6), §§ 387-8 and footnotes. Timoshenko attributes to Vening Meinesz 
(1911) the first suggestion of the soap-film method. 


Hodge also discusses the residual stresses consequent on unloading from the ‘fully 
plastic state’. 
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nvestigation—not yet published—of plastic straining in some systems of 
plane stress and of plane strain. For this a physical hypothesis was needed, 
.dditional to the stress-strain diagram of Fig. 1. The chosen hypothesis 
summarized in § 6) differs in some respects from that of Hodge; but in 
relation to torsion it is equivalent, being simplified by the fact that here 
the directions of shear-stress do not alter in the enclaves (regions of plastic 
strain) as these extend. 

Restriction to the ‘fully plastic state’ is obviated by a recourse to 
‘elaxational’ methods of computation (§ 2). The elastic-plastic interface 
rboundary, which Hodge regards as ‘exceedingly difficult to obtain’, can be 
thus found with an accuracy sufficient to define the plastic stresses, strains, 
ind displacements quite as closely as they need to be known in practice ;* 
nd knowledge restricted to the ‘fully plastic state’ seems insufficient— 
more especially when the section is ‘hollow’ as in the example treated here. 

4. This relates to a solution obtained by Shaw (10). In Fig. 2a small 
numerals give values of his & (= 1500"/D?+ when the sides of the outer 
rectangular boundary are 4-8D and 9-6D) and contours exhibit the gradient 


j 


ff J—except in the enclaves where that gradient has a constant value 


1900 (very nearly). This implies a limiting elastic stress fj; and strain yy 


lated with the twist per unit length (7) and with the dimensions (D) by 
Vy Try 49 
TD prD 15° 
Contours within the enclaves show the distribution of (A—1), a ‘non- 
dimensional’ measure of the plastic strain. (The total (elastic plus plastic) 
strain is measured by Ay,-. Its plastic part (A—1)yy is some sort of 
measure of the ‘work-hardening’ to be expected in a real material.) 

In Fig. 26 the corresponding distribution of the axial displacement w is 
exhibited by contours of the non-dimensional quantity w/D?7. Within the 
enclaves these have been derived from (31), § 13—which is equivalent to 
Hodge’s equation (18). 

Figs. 3, a and b, are corresponding diagrams for a case in which plastic 
straining has proceeded further. The twist (7) is increased in the ratio 6: 5, 


SO now 


3 ( 
yy _ fr _ 49 
7D prD Ls 
and the non-dimensional stress-function ys (of which the limiting gradient 


is 4.900 as before) now stands for 1800 ’/D?. 


SI hi lution which was the basis of Figs. 2 suggest that he judged the 

ndaries of | mputed encla to be known within one-tenth of a mesh-side. I would 

egard this I reasonable, having regard to the fact that ‘Y has continuity, at an 
, botl ids its value and its gradients. 


theory, is defined in (3), § 5. 
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Basic theory 

5. Adopting the ‘semi-inverse’ method of Saint-Venant, at the outset 
we make assumptions which are both plausible and convenient, and we 
accept any restrictions that these entail. The first is that all of the six 
stress-components are zero except X., Y., and that these are independent 
of z (when Oz is the axis of the twisted prism). By it two of the stress. 
equations of equilibrium are made identities, and the third is reduced to 


-X,+ : Yr. 0. 9 


oT Z {2 


X, PT ° r. wet caer oe (3 


In (3), » as usual denotes the modulus of rigidity, 7 the twist of the shaf 
per unit length, and ¥ a ‘stress-function’ still to be determined. W. 
assume both t and‘¥’ to be independent of z, in plastic as in elastic straining, 
To (3) must be added 

X, = Y, = Z, = X, = 0 (4 
(which express the first assumption) and a condition that must be satisfied 
at the cylindrical surface because (6, § 129) ‘shear stress cannot cross an 
unloaded boundary’. For on that account 


X.cos(x,v)+Y_sin(#,v) = 0, 
so pr—=0, ie. ¥ const. (5 


according to (3). For ‘solid’ sections the constant may be made zero, since 
it does not affect the stresses: for ‘hollow’ sections it may be made zero 
on the external boundary, but at an internal boundary it will have a finite 
value, calculable by methods outlined in (13), §§ 189-90. 


The physical hypothesis 

6. Some hypothesis in regard to plastic strain is needed, and (ef. §3 
the standpoint of an earlier investigation has been maintained. This made 
the following assumptions: 

(i) Occurrence of plastic straining is determined by the Mises—Hencky 
criterion. That is to say, normal (elastic) relations hold until the principal 
stresses P,, Po, Pz attain values such that 

(P2—P3)°+ (P3—Py)°+ (Py — Po)” = 8k*, (6 
k being constant for any definite material; and thereafter (in a region 0! 
plastic strain) the relation (6) still holds. (Work-hardening is not contem- 
plated.) 
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(ii) Where the material has yielded plastically the strain has both an elastic 
and a plastic part. The elastic part is related with the stresses by the 
usual equations 


’ 


| P1—O(PotPs)j,--. ete., ( 


~I 
— 


(€1) x E 
in which # (Young’s modulus) and o (Poisson’s ratio) have their normal 


values. The plastic part 1S equivoluminal —that is to say, 
(€;)p+(€2)p+(€3)p 0 (3) 


and in it the principal incremental shear-strains are directly proportional 
to the principal shear-stresses and occur on the same planes. That is to say, 
if 41. Gos Y3 denote principal shear-stresses and Ay,, Ay, Ay; denote principal 


incremental plastic shear-strains, then* 
(Ay,)p: (Aye) p: (AyY3)p = 91:92: %- (9) 


7. The meaning attached to incremental strains can be explained by 
reference to motion opposed by solid friction—slipping that has close 
similarity with the ‘internal slip’ which results in plastic strain. A weight 
made to slide on a rough horizontal table moves at any instant in the 


direction of the force then acting on it, so its incremental displacement 
has the direction of the force; but clearly, if that direction is not constant, 
the same will not normally be true of its total displacement measured from 
some fixed point. To postulate a relation between total displacement and 
instantaneous force would be (in effect) to endow the weight with a 
‘memory’ of earlier slips; and similarly, to postulate a relation between 
total strains and instantaneous stresses would be to postulate ‘memory’ 
in over-strained material. It is the ‘incremental’ principal plastic strains 
that must (in isotropic material) have the directions of the principal 
stresses, in cases where these alter as a plastic enclave extends. 

When they do not alter, neither do the directions of principal strain, and 
hence, in (9), total may be substituted for incremental strains. Then an 


equivalent form of (9) is 
¥1/ 491 Y2/42 3/43 3A (say), (10) 


Yrs Yo: Y3 Genoting the fofal plastic shear strains; i.e. the (shear stress) 
(shear strain) relations have the same forms in plastic as in elastic strain- 
ing—the only difference being that A, in (10), can have any positive value. 


(A = 0 when the strain is wholly elastic.) 


* These assumptions in regard to the plastic strain were suggested by Reuss (14) in 


1930 


The motion is assumed to be so slow that accelerations may be disregarded. 
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8. In the torsion problem Ox, Oy, Oz are principal planes, and of the 
principal stresses 


Ps = 9, ra—fi-—% 1 
where q@* = X?+-Y3; J M1) 
so the criterion (6) reduces to 
3q2 < 4k? (12) 


(the inequality relating to conditions of elasticity), and thus is expressed 
by Fig. 1 with fp = 2k/v3. 

Clearly the directions of the principal strains will alter as a consequence 
of plastic straining—otherwise contours of ‘Y would be invariant. But aj 
the plastic—elastic interface these contours (being continuous) will conform 
with the ‘Prandtl roof’, so they are parallel with the external boundary 
not only after but just before the first occurrence of plastic straining, which 
thus entails no change of their direction. Accordingly it is permissible to 
substitute (10) for (9) as the relation between the principal plastic shear- 
strains and shear-stresses. 


9. The elastic parts are similarly related, except that 3A in (10) is 
replaced by 1/u. So the total strains induced by the principal shears (+q) 
have their directions and are measured by 

+q(u*+ 3a), 
the third being zero since X, = Y, = X, = 0. This means that 
pe., = (14 3Ap)X., Hexy (1+ 3Ay)Y., (13) 


—that is, according to (3), 


Ou } Y 
ao = be = (14 Bu), 
Cz Cx Cy 
(14) 
AX 
pa Wease - (1- 3Au)r— ‘ 
Cz cy : Ca 


cu ov cow Cu ov ~ 
= me == 0 (15) 
CX cy Cz cy Cx 


(the first three being expressions for the principal extensions, the last an 
expression for the shear-strain e,,,). Hence, if rigid-body displacements 
are disregarded, 
u TYZ, v TAX, 

l ow Y (16) 


bs 
= (l- 3Au): +y, — = (1-4 3Ayu) ‘ +2, 
oy T OY Ox 


l ow 


T OX 


7 denoting, here as in (3), (5), and (14), the (uniform) twist per unit length. 
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10. In order that these expressions for @w/éx and éw/éy may be com- 
vatible, the condition 


C v} oy cy 


must be satisfied. Where the strain is wholly elastic it is satisfied because 


Ww " 
“(4 a (aT) +2 0 (A = 143Ay) (17) 


4 =1and'¥ is subject to (1). Within an enclave VY (being determined by 
the ‘Prandtl roof’) is a known function of x and y, so (17) is an equation 
which determines A; and knowing A we can deduce the displacements 
ind the strains from (16). 


The gradient of ‘ within an enclave is directed along the normal to the 


boundary and has a constant value fy/ur = y,;/7. In cylindrical co- 
ordinates the form of (17) is 
ee 1 oA OY 
AV?4 + +2=0, 
Or or r- 0o0 ob 
(18) 
. 79 o* l c ] C2 | 
, ore ror r2 E62’ 


ind if C’, the origin of coordinates, is a centre of curvature of the cylindrical 


boundary, then (cf. Fig. 4) at a point on Cr lying within the enclave 


ew ey ew e2y 
Yy/T; Q. -=—Q, = 0. 
ai om ob oG= 


( l 
so (18) reduces to | A 27/yy 0. 
cl / 


positive sign applying when Cr is directed inward relatively to the 


(19) 


rternal boundary, on which (ef. § 5) ‘’ = 0. Elsewhere in the section 


Y> 0, so oF /@ 0 when Cr is directed outward—a case exemplified by 


the corner enclave in Figs. 2 and 3: then the negative sign must be taken 
19). The intermediate case in which the external boundary is straight 


exemplified by the mid-side enclaves) will receive attention later. 


Computation of the plastic strain 


11. Taking first the case where Cr is directed outward, we insert the 


negative sign in (19) and obtain 
A “(1- rs) — (20) 
a eg ig 
s the solution making A | on the plastic—elastic interface (r = 1r;, say, 
wep 
| ” ; —— C ; 
ere 7, is known). Since the total shear strain is —Ar according to 


cr 


ve it is Ay, according to the first of (19). 


14), at any point in the encla 
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So (A—1) = i: en ae (2) 
r Yy \' 


is a measure of the amount of the plastic straining. 
The case in which Cr has an inward direction is not presented in Figs, 9 
and 3. In it the solution corresponding with (20) and (21) is 


poe oe! ee ee ce | 99 
r Yr \r z 


because the positive sign must be taken in (19). 
Finally when the boundary is straight (r; = ©), by writing 
r=r+é, 
where € is the distance measured in the outward direction from the inter. 
face to a point within the enclave, we deduce both from (21) and (22) that 
A—1l 2ré/yy, (23 
since the positive sign applies in relation to (21), the negative in relatio 
to (22). Thus the plastic part of the shear-strain is now proportional to ¢. 
To summarize this investigation of the plastic strain: It is measured by 
(A—1)yy, which is zero at the interface and within an enclave is given by 
(21), (22), or (23) according as the boundary of the prism viewed from 
outside is convex, concave, or straight. In Figs. 2a and 3a the variation 
of (A—1) is shown by contours: for the corner enclave these were deduced 
from (21), for the outer (mid-side) enclave they were deduced from (23). 


Computation of displacement 

12. Displacements have less practical importance, and it is not easy to 
see how they could be measured (for comparison) in an experiment. (The 
warping of a terminal cross-section could perhaps be measured by an inter- 
ference method, but the theory, since it presumes the cylindrical boundary 
to be stress-free, is not valid near an end.) It is, however, of interest to 
examine, on the basis of (16), the nature of the axial displacement (vw). 

Within the still-elastic region A = 0 and so, by (16), 


l cw O 


u —TYZ, <= +22, = —{¥4} (22+ y?)}, | 
l ow é ‘ ie 
—- = - fy T 4 (x? ty*)) 
7 Cy Cx 
the last two making w single-valued when (1), § 1, is satisfied by ‘VY. But 
‘’ as found by relaxational methods will not satisfy (1) exactly, and hence 


will not permit them both to be satisfied by w. Since they make w plane- | 


harmonic and such that 
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when vy and s have any perpendicular directions, in practice the best 
nrocedure would seem to be to choose some rectilinear closed contour 
yithin the still-elastic region and evaluate w on it by means of (25), dis- 
wributing any ‘closing error’ along its length; then to deduce w elsewhere in 
the elastic region by means of one or other of (24) applied along a ‘string’ 
of the relaxation net. 


In this way w can be evaluated up to and on an interface, and then it 


y 
fairy of Prism 














> 


can be computed within an enclave; for its value is continuous at an inter- 
face, and in an enclave the resultant shear-strain is Ay; and is directed 
perpendicularly to Cr. 

13. Referring to Fig. 4, let the coordinates of C (§ 10) be x, y,, and let 
Cr be inclined at an angle « to Ox. Then the coordinates of a point 7 on 
Cr ; : 

r are x Y~+?r COS a, ¥Y = Yo +rsina, (26) 
and the component strains at the point (if within the enclave) are 

e.. 1+), sin a, iss Ay, cos X, (27) 
according to what has just been stated of their resultant. So by (14) and 


16) of § 9 (in which 1 dAu A) 


ty¥—Ayy sina “ey — 3. 1+), COS a—Ta (28) 
Cy C2 
nd by use of the relations 
( ¢ sIna < c COSa ¢ ° c 
COS a , + sin « (29) 
o7 } Cx cy rf ca cr 
ow cA 
these can be shown to yield the same expression for ——, because — = 0 


Oxucy Cx 
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and (Cr being now directed outwards) 


oA A 
—+-——2r/yy = 0, according to (19). 
cr 7 
Also, since (29) imply that 
cd c ° C 
COS x +sma—, 
er Cu cy 


it follows from (28) that at a point on Cr within the enclave 


1 ew : ; 
-= ycosa—a sina = const. = y, cos a—2,sin x (30 

T or 
if Cr cuts the interface at x;, y;, where w = wy, (say). So at a point on (; 


within the enclave ‘ ; 
w= w,+7&(y, cos a—2,siNn a) (31) 


when € has the same significance as in $11. This last relation can also be 
established in the other two cases (i) where Cr is directed inwards, (ii) where 
the boundary is straight. Hodge’s relation (18) is equivalent. 


Conclusion 


14. In itself the plastic torsion of non-circular cylinders has not much 
practical importance. But Figs. 2a and 3a (which exhibit computed 
strains) perhaps have wider significance for their bearing on the general 
question of ‘relief of stress’. The rather small radius at the interior corners 
would make these regions of rather intense shear stress if the material 
had remained elastic; and their small extent as revealed in Shaw’s deter- 
mination of plastic stresses might seem to imply that these will be accom- 
panied by large plastic strains. In fact, even in the circumstances of 
Figs. 3 (which § 4 shows to entail very considerable twist), the plastic part 
of the shear-strain only just attains a value twice as large as the limiting 
elastic strain yy. Thus in this example only moderate strains are required 
to relieve a quite intense stress-concentration; and no feature of the 
torsion problem suggests that the same will not be true of more complex 
stress-systems. 

For the diagrams which illustrate this paper, and for the computations 
which are presented in Figs. 2 and 3, grateful acknowledgement is due 
to Miss G. Vaisey and to the Clothworkers’ Company, whose grant made 
her help available. The success with which she has treated this ‘hollow 
section gives grounds for believing that any shape can be similarly investi- 
gated (i.e. plastic stresses, strains, and displacements all determined) with- 
out recourse to experiment. 
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SUMMARY 

This problem was analysed some years ago by one of the present authors, using 
strain-energy methods, but the results were too complicated for application. In th 
present paper the spokes are replaced by an idealized disk, a method which has given 
useful results in other problems, e.g. that of a cantilever composed of parallel beams 
interconnected by cross-bars, the critical load of a battened strut, ete. 

The idealized disk is assumed only to transmit actions similar to those transmitted 
by the actual spokes, i.e. it applies bending moments, radial forces, and tangential 
forces to the rim, but neither radial shearing stresses nor circumferential stresses 
occur in the disk. The analysis of the idealized wheel leads to explicit solutions for 
the radial and tangential displacements of the rim, from which the intensities of all 
actions between the disk and rim are calculated. Integration of these actions over 
the area of disk representing a particular spoke is assumed to give the corresponding 
actions for that spoke. 

Examination of actual wheels indicates practical limits for the value of certain 
parameters which permits considerable simplification of the equations. Direct and 
easily applied formulae are thus obtained for the axial loads, bending moments, and 
shearing forces in the spokes, and for the bending moments in the rim. 

Experiments on wheels cut from mild steel plate with twelve, six, and three spokes 
respectively showed very good agreement with calculation from these formulae. 


1. IN an earlier paper (1) strain-energy methods of analysis were applied 


to the problem of stress distribution in a wheel having spokes capable of 


transmitting bending actions. The results, however, were complicated and 
the degree of accuracy required in computation made them impracticable. 
In the present paper the same problem is approached by a method which 
has given useful approximate solutions in other problems (2) and much 
simpler results are obtained. Fig. 1 shows an artillery wheel in which 

R is the radius to the neutral axis of the rim, 

R, is the radius of the hub, 

Ll is R—R,, 

N is the number of spokes, 

A, is the cross-sectional area of each spoke, 
EI, is the flexural rigidity of each spoke, 
A is the cross-sectional area of the rim, 
EI is the flexural rigidity of the rim. 


(Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 4 (1949)] 
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It will be assumed that the spokes are replaced by an idealized disk 
hich can only transmit actions similar to those transmitted by the 


spokes; i.e. the disk can apply bending moments as well as radial and 


tangential forces to the rim, but neither radial shear nor circumferential 
tresses will be transmitted in the disk itself. At radius R, (i.e. at a point 


p Radial 


P Tangential 






































distant R.. from the centre) the disk will have a thickness A, = N.A,/27R, 
and a flexural rigidity per unit length EI’, = NEI,/27R,. The values of 
these quantities at radius R are A and EI’ respectively. The stresses in 
this idealized wheel are determined, and it is then assumed that the 


resultant forces in a spoke of the real wheel are the same as the integrated 


forces in a disk segment centred about the spoke position having an 
ngle equal to that between actual spokes. 
\t an angular distance y measured from any origin, as shown in Fig. 2, 


] 


let the disk apply a couple of intensity m per unit length of rim, a tan- 


gential force of intensity p, and a radial tension of intensity @. 
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The equations of equilibrium of an element subtending an angle 8; at 
the centre of the wheel are, then, 


zo = R(m—F), (1) 
oF 

H+#tR, 9 
Cus “ 
en a (3) 
Cus 


where M, H, and F are the bending moment, tangential force, and radia] 
force in the rim at the section. 
p 








ate 


If uw be the radial displacement of the element measured outwards and 


Fic. 2. 


v the tangential displacement of the element measured anti-clockwise, 


H = ao (1 =) (4) 
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The element of the disk may be considered as a cantilever attached to 
the hub, carrying at the rim junction a moment mF dy, a tension ¢R dx, 
und a shearing force pk do. The rim junction is displaced w radially, 

1 cu v 


Reb R 


tangentially, and rotated by the amount If x is measured 


radially inwards from the rim 





— | 
El oa px—m 
; 6HI'| Cu 
: vhich g eS ) / fi 2 = . 6 
dial | hich giv 7 RE \ ab (2R Nui, (6) 
2EI'(.,éu 
and n 2] L (3R—2))v}; 7 
Ri? | edb sos ei 4) 
: AEu (8) 
i : j . 
If the ring is assumed to be inextensible, from (4), 
ov 
’ ne ; 9) 
' Cus \ 
Equations (1) to (9) give 
wire h e *U ! PA dn Lhe u 0, (10) 
Calf L A /;A 2 Ou? 3 
where ae 
2R?] 
k, = 21 
iT 
$4 Re’ t4A 
k, l 3R—2I)4 
PT am: 
: t{Re]’ ‘ ‘s 
nd h (3R?2—3 RI+/?). 
; 3] 
The solution of (10) is 
i] e ( osh us C’ eos ¥ cosh Bus t Dsin xls sinh Bus T 
B sinh ys S sin ous cosh By | F cos vif sinh By, (11) 


where J, C, D, B, S, and F are constants and a, f, and y are defined as 


follow Ss 





Let 
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Then 


and p= (G47 +34 3 v2 > X fy \> 
-_—” % 9! | 


2. Case 1. Radial load on rim 
If a load P acts radially inwards on the rim at 4 = 7z, wu will be the same 
for equal positive and negative values of y%. The constants B, S, and F 
are therefore all zero and 
uJ cosh ys | C'eos xs cosh Bis |. J) sin ves sinh By; (12) 


0 
- ga («sin vs cosh Bys+-B cos aes sinh Bis) 


Os 
v= ——sinh yb — 
‘ + 


ee ; 
—_— pa (Psin xis cosh Bys—a Cos ox sinh Bys), (13) 
v2 B2 
the constant of integration being zero since v = 0 when % = 0. When 
ub 7, Uv 0, and éu/és = 0; also 
tReos¢ db— pRsinds dibs LP. 
0 0 


These conditions lead to the following values for the constants J, C, D. 


J $P cosech ya[y(y? +-1){(a?-+- B?)?— 2(a?— B?)-4 1} | : 


~ 
_ 
~~ 
S 
— 
= 
to 
— 
& 
tw 
— 
to 
bo 
— 
to 
-. 
-] 
ro 
— 
Se) 
tw 
— 


| 
j 


C : J sinh yz > 


B(a?+ B y?)sin az cosh Ba+a(a?4 p?- y")cos v7 sinh Bx 
vBy(cosh 28a — cos 2a7) 


D I sinh yz : 


B(a?-+-B?—y?)cos az sinh Ba —a(a?+ B?+-y?)sin az cosh Br] 








xBy(cosh 282 — cos 2a7) 
RAE 6EI’ 6(2R—1)EI’ 
where a es; b=—_- 3 FP on E., 
l I? is 
If J’ = 0 the disk has no flexural rigidity and the problem reduces to 
that of the wire wheel with radial spokes. Then 


k=2, k=K+1, 
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¢ and Y = 2 


1f 


vi3{(K+1)?—]}]; 


and the constants J, C, and D deduced from 


them, are those previously obtained for the radially spoked wire wheel (3). 


Ifa tangential load P acts clockwise on the rim at 


Tangential load on rim 


7, radial displace- 


ments are equal but opposite in sign for equal positive and negative values 


and 02a /éds? 


, and D are zero, 


S sin as cosh Bib 


{8 sin ad sinh Bub 


fx sin ays sinh Bys+-B cos os cosh By} |+-G. (15) 


0: also 
us dub 
0 


Rdb—2 | md 
0 
x 1)? B2\5 ( X 
I) (a2 + B2+-y?)? 


.-—y*)COS v7 sinh Ba 


) )/ 
.8(cosh 28a 


vB(cosh 2Ba 


where 
ind 
from which 
M L(A TE V | 
Then 
A l { i 
0, xv = /4{(A 
Sal 
— also 
1d 
The constants a, B, and y, 
3, Case 2. 
13 me 
The constants J, ¢ 
Vhi l B sinh ys 
nd 
[B S 
cosh yi . 
_ eae ) 
F 
D q 2 IZ 
X 2 
When ub 7, U 0 
| tRsin 
| 0 
) - 
und 2 |p 
The constants are then 
P (y*+1 
B : 
2 2 sinh yx (a+b 
S Bsinh yx| : 
F Bsinh yn| ; 
eS 
0 r 17 2+ (a* 


—-| aya l(a h 





F cos ovf sinh Bus 


[ p Reosus dis 


] 


\2 —y*)sin az cosh Ba 


2(3.R? 


(14) 


x cos ax cosh Bis} + 


1P 


e. 





y°B* 
248 sin az cosh Bx 


cos 2a77) 


|- 28 cos am sinh Bx 
cos 277) ' 


3 Rl | 
3 R1+-1?)b | 
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4. Application of results 

Having evaluated the constants the displacements of the rim are com- 
pletely defined and the actions at any section of it or in the idealized disk 
can be calculated. It is assumed that these rim actions will be nearly 
enough the same in the actual wheel with a finite number of spokes and 
that the actions in any spoke can be found by integrating the corresponding 
forces in the disk over the segment replaced by the spoke. It would be 
more logical to integrate the components parallel to a spoke, but since the 
number of spokes in an actual wheel is 12 or more, the difference is very 
small and the slight gain in accuracy is not worth the loss of simplicity in 
equation (16). The axial load in a spoke at %, from the origin is thus, 
provided that yy lies between 6 and 7—8, 


bo 6 AER ho 6 
7+ = tR di - ; u di. 
tho -0 bo—0 
Using equation (9) this is 
= REA 
Ty, - eee 1 (Vy, +0~ Vbo- 9): (16 


In the particular cases of the spokes at x, = 7 and %, = 0 we have fora 
radial load at 4 = 7 
2REA 2REA 


Te : ] (v,-9); to ] (v9). 


For a tangential load at 4 = z, from considerations of skew symmetry, 
7 ‘Zbl 
T, = TN = 9. 


The bending moments in the rim can be found from equation (5) and the 
shearing forces and bending moments in the spokes at their junction with 
the rim from equation (6) and (7) either by direct or graphical integration. 

If the full equations are used for the determination of the resultant 
actions, the integrations lead to lengthy expressions and since in most 
practical cases they are unnecessary, as will be shown in the next para- 
graph, they are not quoted. 


5. Simplification of results 

The formulae deduced are rather elaborate and can with little loss of 
accuracy be considerably simplified by making certain reasonable assump- 
tions. 


Table 1 gives the leading details of a number of artillery wheels and 
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Table 2 the various constants calculated from them. It will be noticed 
that « and £ are very nearly equal and of such magnitude that sinh az 
and coshaz are practically equal; y is so small that sinhyw may be 
replaced by yz, and k, is very large compared with k, and k;. Further, 
) and d are negligible compared ‘with a. 

It will therefore be assumed that k, and k, can be neglected and 


K — R4A/I1. It follows that a = B = (K/4)}. 


TABLE | 


All dimensions in inches 

l 

1g ] A R R 

I 1*455 0°057 0°134 25°50 16°25 0°64 
I 39 2°055 0°260 17°5 11°875 | 0°68 
57 0°130 0°270 29°5 20°25 0°69 

1125 | 0028 0-06 32°0 24°75 0°77 

5 3°76 0°165 0°336 25°0 16°0 0°64 
4°O 0°132 0°275 27°5 18-875 0°69 

12*4 1°40 0°053 Or125 35°5 26°0 0°73 
133 0°00035 0°02955 55 6°5 0°765 


TABLE 2 


Meet e 


N M \ y ' 8 ‘ I ! } 
15 0°017 3°28 3-28 3°29 0° 309 0°0013 0°0026 

6 oO 3°O7 3°04 3°06 0°383 0°0037 0°0072 

r 3°19 2°99 3-08 0°393 o-O019 0°0035 

oO 3°82 3°82 3°83 0°087 0°0003 0°0005 

15‘OI o'l3 2°94 2°73 2°33 0°525 0°0039 0°0053 

75 2°9 2°79 2°55 0°405 0°0022 0°0042 

) oO 4°0 re) 4°0 0-170 0°0005 0°0009 

3°53 3°73 79 0°039 00-0001 


Radial load on rim. With these assumptions the constants for this case 


J i a " 
27 REA 
0 Plae-*" (sin am -+-cos a7) 
REA , 
D Plae-°7(sin am—cos xT) 
REA 
Then 


_ Pl ' 
 2REA\a 


cosh X7T 


/ 


/ cos oas(cos am +-sin a7) —sinh oa sin oaf(cos aa —sin v7)}]. 


{ ] 
COSN a 


When ays is 3, the difference between sinh a and cosh a is only 0-5 per 
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cent.; this difference rapidly diminishes as oz increases. Since « for prac. 
tical wheels is about 3 or more, very little error is introduced by putting 
sinh a equal to cosh op for values of % greater than 47. At this value. 
however, the periodic term in the expression for wu has become so small 
in comparison with 1/z that it may be neglected, and for all values of j 
the expression for u can be reduced to 


—_— sewal- —ae-*?(sin ab+ cos ad) . (17 
where ¢ a—w. Then 
_—_ Pl if | on 
; sama (-*+ wes 4) 8) 
and a= xe x¢(sin uh + cos 4): (19) 


The maximum compressive force in a spoke occurs under the load, and 
if 26 is the angular spacing of spokes its value is, approximately, 
i$ 
T,=2 | tRdy -P(1—N-1—e-“9 cos a6). (20) 
a—0 
The maximum bending in the rim also occurs under the load and is 
obtained from equation (5) when % = 7. Its value when « is large is 
approximately 


Pl 08 
M, — = ‘ (21) 
ReA 
The maximum tensile force in a spoke occurs when % = 0 and is 
8 
TM) = 2| tRdfb = PIN. (22) 
0 
The force in any other spoke at ys = xu, is, from (16), 
T',,, P| N-1—« x0(cOs ad, COs af sinh af -+-sin ad, sin a4 cosh x9)]. (23) 


The intensity of bending moment at the rim junction is given by equation 
(7) which, using the same argument as for uw, reduces to 
PI’ ob 
m= — (83 R—21l)—+e 


RIA xO{4]y2 sin xp — (3 R—21)cos ad}. (24) 


An example of the shape of the curve is shown in Fig. 3 which has been 
plotted for « = 3-79 and / = 0-765R. 

The maximum value occurs very near to the load, and the simplest 
procedure is to plot the curve and find by trial the maximum area on 4 
base length of 26. This will be the approximate value of the maximum 
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bending moment in a spoke at the rim junction. In the same way the 


intensity of tangential force deduced from equation (6) is 









































3PI’ ab a : 
V pa (2Rk—I she, ce xO 9]y2 sin adb—(2R l)cos P| Y (25) 
| “ | 
| P 
| 
‘ | 
| 0 
| | 
| 
| 
5 a 
ay __|-:005 | ce 
| 
Radial oO 
--O10 
ee = am -02 
| 
| | 
O15 
Vl — io | a 03 
i a | | mR 
: | P 
’ - 2 130 10 
¢ Degrees 
Fic. 3. 
This eurve is also plotted in Fig. 3 and the maximum value is very near 


to the position for maximum m. The curve of p can be plotted and the 
maximum area on a base 20, found by trial, is the approximate value 
of the maximum shearing force in a spoke. 
Tangential load on rim. On the same assumptions as before, the con- 
stants for the tangential load reduce to 
B Pl Ple-%7 sin ar : Ple-°* cos am 
3 


27REAy = =———t~—<“<i~ti‘“SCSCSRRiENSC’ ~ PBA 
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Q- Pl {_ 1 _ 12R 
~ QI7E|RAy? 4/1'(3R2?—3RI+P)° 
Then 
2 
2 aes (t—< 2¢ COs t), (26) 
_ PRP PL fee ‘ wy) 
~ Sn ET'(3R2—3RILE) + SREA| a OS? “p)— > (22) 
> 

and t ont x? COS ). (28) 


The load in the spokes at % 0 and & m7 are zero from considerations 
of skew symmetry; the load in any other spoke at % = x, is, from 
equation (16), 


r P| i x 


Yo N 2a 
< {sinh a cos a6(cos ady—sin ady) + cosh af sin «6(cos ady+ sin r$0)}]. (29) 
From equation (7) 
21'P{1 . Pl 3R—2l 
m = — —_ xe~*?(cos ad-+sin vg) | — 5 513 
RA |x }  4n\3.R?—3RI+-/? 
but the first term is very small in comparison with the second and can 
be neglected, so that approximately 


Pl 3R—2l ; 
m= — . (30) 
4 \3.R°—3 RIP? 
Similarly, from equation (6) 
3P 2R—I ; 
p= —- =|. (31) 
47 \3R?2—3RI+(? 


The curves of t, m, and p for the case when «a = 3-79 and / = 0-765R are 
shown in Fig. 4. 

The maximum bending moments and shearing forces in the spokes at 
the rim junction are 


IP € . 
vw’ — —s | 3R — 2] (32) 
2N \3R?—3RI+0F* 
3PR 2R—l 
‘ Fr — — (33) 
wn 2N (; R2—3Rl TY 
The bending moment in the rim is, from (5), 
M = spn e $(cos ad+ 2x? sin ~6)]. (34) 
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The maximum value of this is best determined from a curve such as shown 


in Fig. 4. The maximum value for this case occurs approximately at 








| Tangential load 
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6. Verification of results 

To test the reliability of the results an experimental wheel, No. 8 in 
Tables 1 and 2, was cut from steel plate 0-355 in. thick. This was loaded 
both radially and tangentially, and the forces in certain spokes were 
measured by a mirror extensometer. The moments in one spoke were 
btained by a M.I.T. moment-indicator to check the calculated moments 
and shearing force. The tangential movements produced by a tangential 
load were also measured. The experimental results are given in Table 3 
together with the calculated values. 


rhese tests are not of equal significance; the tangential load was limited 
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by stress considerations to 400 Ib., which produced a maximum spoke load 
of about 90 lb. with a direct stress of only 660 lb. per square inch. Thy 
accuracy of readings at such low stresses is doubtful and the measured 
values of 7’ were erratic; a spoke next to the loaded point gave exagt 
agreement with the calculated value, while the corresponding one on the 
opposite side gave a value 18 per cent. too high. The numerical value of 
these should be the same, and since the measurements are equally reliabl 
it may be assumed that the actual and calculated values differed by about 
9 per cent. 
TABLE 3 
12-spoke wheel 


Load Iclion Position Calculated 
i measured bo Experiment | Calculated | Error ‘ from equatior 
(| —0°745 0-7 16 
Radial T/P 7 i _ 
o-714 3°5 20 
TIP or 0220 0°221 o 29 
T/P On 0270 0:221 18-0 29 
v/P 7 0:000095 in. | 0:0001 in. 5:0 27 
Tangential 4 v/P ° 0000097 in. | O°0001 in. 30 27 
W’/PR (rim) 7 0°0405 0'0365 10'0 32 
} } i hs 
| , 
| | M’/PR (hub) 7 0°055 0°055 ° 32 & 33 
| 
F’/P 7 0128 0*120 6:0 22 


The safe radial load was much larger and the measurements corre- 
spondingly more reliable. The maximum spoke load agreed well with 
values calculated either by the general formula of equation (16) or the 
approximation of equation (20). 

Tangential movements could be obtained very accurately by mirrors 
and agreement between the measured and calculated values is close. 

The bending moments and shearing force in the spoke as measured bj 
the moment indicator agreed well with the calculated values. 

A further check on the results was made by Mr. W. C. Brown, who 
analysed wheel No. 4 in Tables 1 and 2 by a relaxation process. 

Under a vertical load the maximum force in a spoke was calculated to 
be —0-73P, i.e. the mean of the values obtained from equations (16 
and (20) respectively. The force in the spoke opposite the load, i.e. at 
% = 0, was estimated as 0-0785P, which is less than 6 per cent. below the 
value of 0-0833P obtained from equation (22). 


The maximum bending moment in the rim from equation (21) 1s 
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Ike Io 0-:0653PR; the value obtained by Brown is 0-0715PR, a difference of 
th. T under 9 per cent. 
easury On the completion of these experiments alternate spokes were cut and 
e exact | after further tests alternate remaining spokes were cut. The results of 
on t experiments on these six- and three-spoke wheelsare given in Table 4 with 
alu calculated values of the forces, etc., for comparison. 
relia ii 
TABLE 4 
Vy and 
6- and 3-spoke wheels 
Val té . 9 , , 
, Caiculaled 
Experiment | Calculated | Error ‘ from equation 
TT 857 852 fe) 20 
c 205 0°000201 20 27 
, 1 82 73 I1‘O 32 
T Its III 3°5 32 KX 33 
= 26 0°24 8-0 32 
| 7 oo! 0°945 4°0 20 
036 0°00040 II‘o 27 
7 135 0146 5°5 2 
7 1of 15d 4°I 32 X 33 
T 436 45 10°00 3.3 
The close agreement between measured and calculated values in wheels 
with so few spokes was unexpected and indicates that the formulae 
col obtained in the paper may be used with considerable confidence for all 
) wit wheels of practical dimensions. 
§)i 
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ON THE NATURAL FREQUENCIES OF 
SUSPENSION CHAINS 
3y A. G. PUGSLEY (University of Bristol) 


[Received 14 December 1948] 


SUMMARY 

The equations for the small oscillations of a uniform suspension chain are intract- 

able. The writer, seeking direct expressions for the natural frequencies of such a 
chain and guided by the results of some elementary oscillation experiments, presents 
a simple approximate theory of the oscillations. Its results are compared with thos: 
deduced by Routh for a cycloidal chain. Some semi-empirical formulae for the first 
three natural frequencies are then put forward and comparison is made with 
experiment. 
1. THE Severn and other suspension bridge projects recently under con- 
sideration in this country have led to a revival of interest in the theory 
of such bridges. Moreover, the failure under oscillatory conditions of the 
Tacoma Bridge, U.S.A., has emphasized the importance of an under- 
standing of their vibration properties. These considerations have led the 
writer to a study of the stiffness and vibration properties of suspension 
bridges, as a preliminary to which some elementary experiments on the 
relevant properties of uniform chains and strings were made. Of interest 
in the present connexion are some experiments on the natural frequencies 
of a chain when suspended between two fixed points on the same level and 
oscillating in its own plane. 

The reasons for making these particular experiments arose out of a brief 
study of the theoretical position as stated by Routh (1).+ Let x, y be the 
rectangular coordinates of any point P on a suspension chain when hanging 
in equilibrium, x being measured horizontally; and let s be the distance 
along the are of the chain to P from some convenient origin, and m the 
mass per unit length of the chain at P. When the chain is executing a 
small oscillation suppose P, at time ¢, to have coordinates «+, y+7; 
then Routh states the equations of motion in the general form 


i. 1 d pe ype (1) 
dt? m ds ds ° ds 

d*» - ] d . dn aa dy (2) 
dt? m ds ds ° ds]}’ 


+ He refers in a footnote to some earlier work by Rohrs (2) in 1851. This has been 
examined, and treats of the equations of motion for the symmetrical modes of a chain that 
is nearly horizontal; results are obtained for the first two modes and can be shown to 
agree with Routh’s when taken to the limit corresponding to d/s —> 0 in equations (16) and 
(18) of this paper. 
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where 7' is the equilibrium tension in the chain at P, and U is the increase 
n tension there at time ¢ when the chain is in motion. If the chain is 
inextensible, its geometry at any instant must conform to a third equation 


da dé dy dy ; 


ae oe 0. 3 
ds ds ds ds (3) 


These three equations, which provide a basis for the determination of &, 
.and U in terms of s and /, are intractable for the uniform chain with 
m constant, and Routh adopts a condition more amenable to solution— 
that of a chain with a mass distribution varying along its length so that 
the chain hangs in the form of a cycloid. In this case solutions are possible, 
though they can only be stated in convenient explicit form when the chain 
hangs with a dip that is small compared with its span. Two related needs 
thus appeared—one for a simple approximate theory and the other for 
explicit expressions, if necessary empirical, for the natural frequencies of 
suspension chains. The aim of this paper is to contribute to these needs. 

2. The oscillation experiments were conducted with uniform strings, 
vires, and chains with various masses per unit length and with spans up 
toabout 10 ft. The ratio of span to central dip was varied from about 5 
to 15. Excitation was either by a D.C. motor via an eccentric and weak 
spring or by hand, and the frequencies, which were recorded by the 
counting of oscillations, were observed to an accuracy of within +3 per 
cent. Only the first three natural modes were explored. The conclusions 


lrawn from these experiments may be summarized thus: 


1) the natural frequencies were independent of the mass per unit length 


of chain; 


2) the natural frequencies were, within the accuracy of the observations, 
practically independent of the span of the chain; 

(3) the natural frequencies varied approximately inversely as the square 
root of the dip of the chain, but in each case were rather more 
sensitive to dip than indicated by the square root law; 

(4) the ratios of the natural frequencies in successive modes were 


roughly 1-0, 1-5, and 2-0. 
The strings or chains being effectively inextensible, the first three modes 


corresponded to small oscillations with one, two, or three nodes in the 


length of the chain, not counting the fixed end points. 


Since this paper was written, my attention has been drawn to a recent paper by 
Reissner (3) in which the Rayleigh—Ritz method is used to obtain an approximate solution 
r the lowest frequency of a suspension bridge with deck girders of low stiffness. The 
inalysis, which is necessarily rather lengthy, leads incidentally to formulae from which 
Reissner deduces ‘the remarkable general result’ that the lowest frequency obeys the 
mditions (1), (2), and (3) here 
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3. A natural way of coordinating the foregoing results was to cast them 
in a form suggested by dimensional theory. The parameters involved 
comprise the span JL (or the chain length s), the central dip d, the mass 
per unit length m of the chain, and the acceleration due to gravity g. We 
may thus write for a natural frequency n that 


n = $(L,d,m,q), 


and thence by the usual dimensional arguments come to the form 


q\+_,(L 
= (a8 


Here the absence of m and the presence of (g/d)! conform approximately 
to our experimental findings, which also suggest that the function  con- 
cerned makes L/d comparatively ineffectual. This ratio could equally well 
be replaced by (s/d), a form which corresponds to the practical operation 
of changing the dip and span without changing the chain length—the way 
in which LZ and d were commonly varied in the experiments. 

4. Another way of approaching the problem is to regard the oscillations 
of the chain as arising from the propagation of transverse waves along its 
length. Then if the velocity of propagation v is treated as constant, the 
natural frequencies will be given by 

UV 
n=1-, (5) 

s 
where s is the length of the chain and 7 has values of 1, 1-5, 2 for the 
successive natural modes. Now the velocity v in a chain under tension T 


is given by T 
v . (6) 
m 


and since 7’ does not vary greatly along a catenary with span to dip ratios 
in the practical range (say 5 to 15), we may expect (5) to give us a fair 
approximation to the desired frequencies. 
The tension at any point in a catenary is given by 

T = mgy, 
where y is the ordinate to the curve measured from the directrix. At the 
lowest point of the chain, y = c, the parameter of the catenary, and to 
a close approximation the average value of 7’ acting in the chain will 
therefore be 


T' = mg(c+-4d), (7) 
where d is the dip of the chain. Now the parameter of a catenary is given 
by (4) 92 

c= —-— hd, 


a 
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them ind inserting this in (7) gives: 
rolved = mgs” || 4d? 
olved ij YJ [)-; aA, (8) 
. Mass Rd N { 38° 
y. We —_ ,; 4d? . ‘ 
Substituting from (8) in (6) and noting that 3,3 8 small compared with 
os* 
unity, 1 fa\4 a2 
; 1 22 
v (2 =s(1——], (9) 
2V2\d 387 
: jas D]2 
~ t q\s at 
whence trom (5) n (! a =}. (10) 
4 2v2\d 38° 
— This expression is at least approximately consistent with our experiments 
ne in that » is almost independent of s (and L), that » varies nearly as the 
yw Con = 
— reciprocal of .d, and that the frequencies of successive modes are in the 
— ratios of 1, 1-5, 2,.... 
ratlor ; : . 
ee We have so far assumed that v is a constant along the length of the 
chain. In fact at the centre of the chain, and to a close approximation 
; for a large central part of the length, 
itions 
Ng ts | q\4 2d? 
» (7 oe (11) 
t, the ZNZ d, q - 
The velocity near the ends rises and at the ends becomes 
5) l /q\4 4d?\ 
, J)? (14. ), (12) 
: 2V2\d ' 
r th / j 
ion 7 | These small variations will probably have little effect upon the frequency 
for high modes with wave-lengths short compared with the chain length, 
6 but for the first two modes we may expect the low value of v over a large 
central part of the chain to increase the coefficient of d?/s? in (10). 
ratios 5. It is instructive to examine the probable variation of the coefficient 
a lal of d*/s* in (10) by reference to Routh’s solution for a cycloidal chain. 
This can be rendered explicit when d/s is small. For the lowest natural 
frequency, Routh gives ’ 
rs Any = 7, (13) 
where A = 2,/{1+-(bk?/g)} and sina, = 1/46. Here b is the radius of the 
t the 


generating circle for the cycloid, / is the semi-length of the chain, and & is 


d to 4] : " . 
- the circular frequency required. For the case when ap, the slope of the 


. will chain at its ends, is small, (13) can be written: 

7 l gq [427b* 
: Ny 4 ( ee a | (14) 
riven 277 A \b\ /? | 


In this expression we may note, in terms of our notation in (10), that 


+s, and to a close approximation, if d/s is small, b = s?/32d. 
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Substituting these results in (14) and neglecting terms in d/s involving 
higher powers than the second, we have 


1 /q\4 32d? 
Ny 9,/9 (“)* (1 — a (15 
2V2\d (77s) 


The frequencies given by Routh for successive modes, cast in the same 
form as (15), are: 


1-43/9\3 /, 32d2 , 
ae sot 3)> 6) 
“ 2v2\d (1-4327)2s2]’ ( 
2 /(g\3 32d? : 
"3 = sald) \ ~ Qa3s2)’ (17 
9.5 1 99,72 
and n oo 19)\8 {5 __ 32d io 
. 2v2\d (2-57r)2s? 


Thus for the higher modes the coefficients of d?/s? decrease and for ng, for 
example, approximate to that in our equation (10). These results, based 
upon Routh’s solution, of course, apply strictly only to a chain with a 
mass distribution such that the curve adopted by the chain is a cycloid. 
However, for a flat cycloidal are when d/s is small no large departure from 
the conditions appropriate to a uniform chain of similar small dip would 
be expected. 

6. The solution for the cycloidal form thus confirms our general result (10 
forthe case whend/sissmall, but some theoretical confirmation of the general 
form of (10) when d/s is rather larger seemed desirable. As a possible guide 
on this point, the case of two equal concentrated masses placed at equal 
intervals along a weightless string of length 3a fixed at its ends to two points 
less than 3a apart was examined. In this case, if the string is inextensible, 
only one mode is possible; each mass moves at right angles to the direction 
of the outer portion of string attached to it, the two masses at any time 
moving in opposite senses. The natural frequency of this system is found 


n seal (4) se | 2 cos), (19) 


where d is the ‘dip’, i.e. the depth from the horizontal line joining the ends 


to be given by 


of the string to the level adopted by the masses, and « is the slope adopted 
by the outer portions of string. 


The third term here, ,/(1-+-2 cos*«), decreases continuously with increase 


9,72 
— P Pa-\ . — 
of dip, as does the corresponding term {1——-] in (10). In the limit, 
3.L? 
when a = 37 and d = a, the two masses hang on vertical strings and 


‘ 1 /g 
Ina d 


as for a simple pendulum. 
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To compare (19) with (10), we take « as small and note that a = d/a 


and 5 od. 


3/q\i 9d2 
Then (19) becomes n 5 (9)(0- ss) (20) 


When the rates of change with d?/s? given by (10), (15), and (20) are com- 
pared, it is clear that our two-mass example confirms the indication of (15) 


based on Routh’s cycloidal solution that the influence of d?/s? is greater 
than indicated in (10). 
7. Having regard to the foregoing results, and to the approximations 


involved in each case, it appears reasonable to adopt as an expression for 


the frequency of the first mode 


l /g\3 3d? 
n el he ee, 21 
; 2v2 (‘) ( s? ) ony 


Let us examine this in relation to some of the experimental results. The 
table given below gives a series of frequency observationsy for a chain of 
given length ; the results for the first mode only are given. 








Frequenc y (m4) 


(in.) Calculated | Observed 

54 1°82 "79 

54 2:08 2°09 
2°20 2°29 


54 
The calculated values are based upon equation (21) and show good 
reement with the observed frequencies. 


y 
1 


On the evidence provided by (10), (16), and (17), we adopt similarly 


», ~ 24 (3) ( -~F}. (22) 
; 22\d; i 


9 a\1 70 2 
ind for the third: No - (3)° (1 mit =} (23) 
2v2\d) \ a" 


Comparison is again made with experimentt in the following table. 


for the second mode 





I requencle 
Ng Ns 





culated Observed Calculated Observed 








2°76 2°65 3°99 3°96 

2°38 2°40 3°46 3°33 

2°10 2°12 3°09 2°9gI 
a ee ; : ; 3d? 

Kesults for relatively low L/d ratios are quoted as a test of the correcting term (1 —- r) 
2] 
t The writer is indebted to two students, P. B. Morice and D. R. Sims, for the observa- 
tions in this table and for the main series of experiments referred to in para. 2 above. 
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Our equations (22) and (23) predict n, and ng satisfactorily, though 
their accuracy is not quite so good as (21) in the case of n,. This may, "7 
however, be apparent only, as the higher observed frequencies, owing to 
the method of counting, are not so accurate as the lower ones. 

Equations (21), (22), and (23), the simple form of which was set by the 
elementary wave theory of § 4, thus appear to provide fairly accurate 
semi-empirical expressions for the natural frequencies of a uniform sus- 
pension chain. 


T 
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THE RESISTANCE OF A SUBMERGED CYLINDER 
IN ACCELERATED MOTION 
By T. H. HAVELOCK (King’s College, Newcastle-on-Tyne) 
Received 8 March 1949] 


SUMMARY 
The problem considered is the resistance to motion of a circular cylinder at a 


nstant depth below the free surface, in particular when the motion starts from 


rest and has uniform acceleration. The resistance is expressed as the sum of two 


terms; one corresponds to the wave resistance for uniform velocity, and the other 


ay be taken as giving an effective inertia coefficient, the variation of which during 
motion is of special interest. The expressions are carried to the second order 


f approximation and have been reduced to forms suitable for numerical computa- 


tion. Curves are given showing the variation of both parts of the resistance during 


e motion, for various values of the acceleration. 


1. For the steady motion of a submerged body with velocity c parallel to 
Ox, the condition at the free surface of the water is 
07076 /0x?+-g 0db/cy 0, 

where ¢ is the velocity potential, Ox is horizontal, and Oy upwards. For 
small values of c this becomes formally equivalent to ¢¢/éy = 0, while for 
large velocities the corresponding limit may be taken as ¢ 0. The same 
effect may be seen if we consider the expressions for, say, a moving point 
source at a given depth below the surface; it is easily seen that in the limit 
the image system becomes a point source for small velocities, while it 
approximates to a sink for large velocities. Some discussion has arisen as 
to the appropriate surface condition to use when estimating the effective 
inertia of submerged or floating bodies; but any argument based on steady 
notion assumes a state which has been uniform for a long time, and 
cannot be applied directly to accelerated motion or motion started at a 
given instant. In a previous paper (1) expressions were given for resistance 
in accelerated motion, but no case has hitherto been worked out. It can 
be seen from equation (30) of that paper that, if we proceed only as far 
as the first approximation, the total resistance separates into two parts, 
the wave resistance and the inertia resistance; further, the latter part is, 
to that approximation, the same as for motion under a free surface 
neglecting gravity and thus corresponding to the surface condition ¢ = 0. 
To obtain a more accurate result it is necessary to proceed further in 
the approximation to the solution. In the present paper we consider the 
problem of the circular cylinder moving at constant depth below the 
surface, examining, in particular, motion with uniform acceleration starting 
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420 





T. H. HAVELOCK 


from rest; the solution is carried to the second order of approximation. 
[t has been found possible in this case to reduce the expressions to forms 
which are not too difficult for numerical computation, and curves have 
been drawn to show the influence of the acceleration upon the resistance 
and upon the effective inertia coefficient. 


2. We shall construct the expressions by the method used in the 
previous paper. With the origin O at a depth f below the free surface. 
Ox horizontal and Oy upwards, suppose a singularity of order n created 
at the origin at time ¢ = 0, maintained for a short time 57 and then 
annihilated. To satisfy the condition at the free surface during this 
impulsive motion, we have for the complex potential function 
9 = Pi2z-"— PR(z—2ef )-, (1) 
where P, may be complex, and P* denotes the conjugate complex 
quantity. 


u 


To obtain the surface elevation in a convenient form we write (1) as 
* 8) 
— 1 n ‘" k i n ; 
Wy) = ( ~) _P KM —leike dye — — (2) p* [ Kr —le—ikz—2Kf die (2) 
‘=. ’ n . “ 

(n—1)! (n—1)! J 
0 0 
The result of the initial vertical velocity acting for a time 47 is to leave 
the free surface with an elevation » given by 


L 





E R 266)" Ps . F -n —ixx—kf dq 
Ui) ve ain Ke ak, (3) 
(n—Il1)! 
0 


where Re denotes the real part. 
The potential function for the subsequent fluid motion due to this 
initial surface elevation, without velocity, is 


Qqhin ae : : 
w— qJ 7 4 oT Ke bp KZ ~2«f sin(gttx?) dk. (4) 
n H 


. 


0 
We now consider this to be a continuous process occurring as the origin 
moves parallel to Ox with a velocity c, with ¢ and P, functions of the time. 
Let s, be the distance travelled by the origin from the starting-point; then 
in (4) we replace t by t—7, and z by z—(s,—s,), so that z is now referred 
to the moving origin. Integrating from the start up to the instant ¢, we 
obtain for the complex potential 
t L 


> I* sN+1,y3 ' 
a a | P&(r) dr [ Le-txe-20r idx; (3) 


w —— —— - 
0 0 


ait 


(z—2if)" (n—1)! . 


with L=e- i{«(s,—s7)-gtxt(t—7)} — e—ifn(s;—87) +94 xht—7)} 
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This result may be confirmed by using the pressure condition at the free 
surface, when the axes are moving parallel to Ox with velocity c and 


acceleration ¢. For these relative coordinates we have 
, ; 
) CO C ‘ a 
. Por —c—— 3q*—gy, (7) 
p ot Ox “, 
lcd ced 


get g cx 


The condition that p is constant on the free surface leads to the condition 


{o-w d*w dw oO [dw 
vu e \ U ¢é — | = 0: - f 9 
| ot? dz? 7 dz 20 (7) y= (9) 


It may be verified by direct substitution that (5) satisfies this condition. 


3. Suppose a circular cylinder, of radius a, centre at the origin, is 
moving horizontally with velocity c. We assume that the potential can 
be expressed as an infinite series of terms like (5) for integral values of 7; 
ind the quantities P, are to be determined from the boundary condition 
m the circle |z a. If we write 

F(x, t) > (—1)"P,(O)«"1/(n—1)!, (10) 
l 


ve have the general expression 


7 P (t)- ( F* «tye ixz—2kf dk 
0 t 
igt | F*(«,7) dr | Le-**-2Sict de. (11) 


0 0 


We may expand the second and third terms in positive powers of z in the 


neighbourhood of the circular boundary, and we get w in the form 
> (P,,27-°+@, 2"), (12) 
are hk 
vith 
Q) | F¥ (ic, there die + 
n! 4 : 
0 
j)r+igt ( _ 1 = ‘ 
—_ F*(«,7) dr Ln" +te-2"F de. (13) 
n 
0 0 


The bound ury condition on the circle gives 


P. ca®+-a?Qt; P. = a™@. (14) 


1 
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Hence, for the quantities P, we have the infinite set of equations 


‘i ; 
P,(t) = ca*—1a? F («, t)ke-**F dx—g*a? F(x,7) dr L*xte-2F dc. 
0 0 Fi 
P(t) = en - F'(, t)e"e-24f dk +- 
n! 
0 


x 
yn : Igtg2n 


ees . F(«,7) dr [ L*K"+4e-2F dic. (15) 
n. } 


« ry 


0 0 


4. We shall only attempt an approximate solution of these equations 
as far as the second order, that is, up to n = 2. It may be noted that the 
condition at the free surface is satisfied exactly, but the condition on the 
circular boundary is only satisfied approximately to the order indicated, 
For the forces (X, Y) on the cylinder we use the general expression suitable 
for axes moving with the cylinder (2), which is in this case 


X—1Y L pt | (ie) |-aparé - pig [ w* dz*. (16) 
dz C 


We shall find it convenient to take the corresponding resistance in two 
paris; thus, to the present order, 


. { (dw\? 4a 
R, Re| “Lei | (“-) de} = —*PRe(P, PB (17) 
az a 
R, —paré- Re| pi 5 [ w* d:*) = —7pa*é+ 2p Re{ = Py]. (18) 
Cc Cc 
Further, from (15), P, and P, are given by, 
t : 
P 2 a® P ia p $2 j f ;P. P. Ld { * cle—2F dc 
_=o — pp \— 573 »—g*a (—tP,(7r)—«P,(r)} dr | L*¥Kie-**T dk 


0 (19) 
t 


. 


x 
4 274 ° 
va 3a ‘ ‘ ae 
er sf3 = lf vals pigtat | (—tP,(7)—«P,(7)} dr | sii 
2 SF: yes 2 2 


0 0 (20) 


If we neglect gravity, we have approximately 
“ f=) * . 


P, = ca?(1—a?/4f?); P, = —ica/sf?. (21) 
From (17) and (18), R, is zero and 
R, = mpa*é(1—a?/2f?), 22) 


the coefficient of 7pa?é being, to this order, the effective inertia coefficient 
for a free surface, neglecting gravity. The next step is to use these first 
approximations for P, and P, in the integrals in (19) and (20) and so 
obtain the next approximation. 
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5. Before proceeding, we may confirm this process by applying the 
method to the case of uniform velocity c, for which the results have 
previously been obtained by direct consideration of steady motion. We 
require the limiting values of the integrals in (19) and (20) for ¢ becoming 
infinite, the quantities P being constant. Putting in the appropriate form 
for L from (6), we have, for instance, the integral 


( dr ( fpitxe—gtxt\t—r)__ pi(ne+gtndyt Micle-2F dye, (23) 


Integrating with respect to + and taking the limiting value of the 
integral in « for t > 0, it is readily found that (23) has the limiting value 
2g' ky c-*[ we-*+- t{a-1—e-*Eii(x)} ], (24) 


where Ky = g/c”, a 2«, f, and Ei is the exponential integral. The similar 
integral with «? in place of «x? converges to 


292s c-*| we-*+-ifa-2 + a-!— e-*Ei(a)} |. (25) 
The integral with a factor x is not required at this stage; being factored 
by P,, it clearly does not enter into the second-order approximation. 
Hence, to this order, (19) and (20) give, for uniform velocity, 


, a* — — . 
P. ca" i P, 21s cat mre . t+-afo-! = *Ei(x)}], 
. oo eon - . 
F. ' f FP, Ky a4 P| xe ¥4-sa-*-+-a*—¢ ~*Ei(«)}}]. (26) 
Rf: ¢ 


The resistance R, is zero in this case; and from (17) and (26) we obtain 
R, Aor? p07 ate—™| 1 — 2x2 a?{a-2 + 2a! — 2e-* Hi(x)} ]. (27) 

This result agrees, to the second order, with the more general expressions 

obtained previously for the wave resistance at uniform velocity (3). 


6. Returning to the general expressions (19) and (20), we shall examine, 
in particular, motion with uniform acceleration y, starting from rest; thus 
we have ¢ yt, 8s = 4yt*. The first approximation to P, is ya*t(1—a?/4f?), 
ind it is sufficient for the next stage to put P, = ya?z in the integrals in 
19) and (20). Hence, to the required order, we have 


t L 
F yart = P,+ig'aty [ tdr [ L* xie-**F de, 
0 0 
t 2 
F : it sgtasy tdr [ L*xie-*F dr, (28) 
S/ , . 


where L* = etlixy*l-7* gtxi(t- 7)}__ pifdeyl—r)+9hnbd—2)}. (29) 
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We now reduce the integrals to a more convenient form. The integration 
with respect to 7 can be expressed in terms of Fresnel integrals ; after some 
reduction we obtain the result 


t 
| L*r dr = p~etv'-a"{q P(pt—q)+-qP(q)— hie} — 


0 
—p~e i(pt +@*Sq P(pt+ q)—qP(q) = bie ~ig*) (30) 
where 2p? = xy, q = g/2y; 
u 
and P(u) = C(u)—iS(u) | et du. 


0 
For the integration with respect to «, we change the variable from « to r, 
given by x = 4gv?/y*t?; and we obtain finally 


t 0) 

cdr L*xie-**F dx = A,+iB, = k*q-! fi A+iB)v2e-88* dy, (31) 
0 0 

t oa o 
cdr | L*Kie-**S dx = A,+iB, k18g-4¢-5 | (A+iB)vte-8he* dy, (32) 
0 0 0 


Oo eh ; "1, 242, ae Ce tas ‘ 
kh? = 2g/y; B = gf/y*t*; Py k(v—3); Po = k(v+4); 
A C(p,)cos p?+ S(p,)sin pi — p3+S(p.)sin p24 
+{C($k)—k-1 sin 1k?}(cos p?—cos p?)+ 
LESak) +-k-1 cos }k*}(sin p?—sin p32); (33) 
B = C(p,)sin p3—S(p,)cos p?+ C(p,)sin p3—S(p,)cos p?— 
—{S($k)+k-1 cos tk?}(cos p?—cos pz) 
+{C($k)—k- sin $h?}(sin p?—sin p2). (34) 
7. For the resistance, we consider first the part R,. This could be obtained 
to the second approximation, but it was thought sufficient meantime to 
examine only the first approximation. The general effect of the second 
approximation is known in the case of uniform velocity; it consists in 
increasing the value somewhat at lower speeds and diminishing it slightly 
at higher speeds. From some rough calculations it appears that the effect 
in the present case would be similar; but for a general idea of the effect of 
acceleration upon R,, which reduces to the wave resistance for uniform 
velocity, it is sufficient to take the first approximation. From (17) and 
(28), we have 


R, = 2npygtattA, = 1287gpa*kB?(a?/f?) [ Avte-8Be* dy, (35) 
0 


in the notation given in (33). 
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For the second part of the resistance we include second-order terms: 


from (18) and (28) we have 


R, 7 pa*y +-27pa*y(1—a? 4f*) -2npygiatc B, ot. (36) 
From (31) and (33), this leads to 
R, 7pa*yp, 
vith p | —(4+-32kB*b)(a?/f?), 
b | (—3- 16Bv?)v? Be-8B* dv. (37) 


Numerical computations have been made for the integrals in (35) and 
37). The quantities A and B depend only upon the acceleration, while 
the instantaneous value of the velocity enters through f. The integrals 
were calculated for two different accelerations, and for about a dozen 
i to 40. For small values of f it 


yas necessary to go as far as v = 4-0 or further, but subdivisions of 0-1 


ilues of B in each case—ranging from 


for v were usually sufficient. For large values of B the necessary range 
for v was less, but subdivisions of 0-02 had to be taken, especially for the 
arger values of k. For various reasons it was difficult to obtain any high 
legree of accuracy in the final results; but it is considered that the 
ilculations are sufficient to show the general character of the effect of 


eceleration upon the resistance. 


8. Some of the results are shown in the curves of resistance. These 
urves show the resistance for a particular value of the ratio of the radius 
f the cylinder to the depth of its centre, namely the value given by 
?/f? = 0-1. We have chosen to graph the curves on a base of velocity c, 
r yt, the abscissae being c/(gf)'. This was partly so as to bring into the 
liagram the wave resistance curve for uniform velocity; this curve is 
shown as R, in the diagram. 

Taking the resistance R, first, the curve A, shows its value for k? = 97/2, 
or for y/g = 0-1418; while the curve B, is for k? 7/2, or for y/g 1-276. 
The effect of greater acceleration is shown in the lower maximum wave 
resistance and the higher velocity at which it occurs compared with the 
urve R, for uniform velocity. It should be noted that if we had graphed 


u 
the curves on a time base, the abscissae for curve B, would be reduced 
to one-ninth compared with those for A,. 
rreater interest. In 


We turn now to the resistance R,, which is of g 


general, the relative magnitudes of R, and R, depend upon the two ratios 
yganda/f. In the diagram, the curve A, shows the resistance R, for the 
case y/g = 0-1418, and a?2/f? 0-1: the total resistance in that case is 
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given by A,+A,. It is seen, from (35) and (37), that the part of the total 

resistance which is simply proportional to the acceleration is 
mpay(1—a?/2f?), 

If we define the effective mass as the coefficient of y in this term, then 

the inertia coefficient is the same as for a free surface neglecting gravity, 

We could, on the other hand, divide the total resistance by y and so define 


T T T I T T T T T T |! an 
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the effective mass at each instant. However, from the way in which R, 
and R&, arise and from their variation, it seems convenient to refer to R, 
as the wave resistance and to regard R, as the product of the effective 
mass and the acceleration; with this convention, the inertia coefficient for 
motion with uniform acceleration from rest is given by the quantity p 
of (37). It can be seen from (37) that p converges to 1—a?/2f? for both 
c—+>0O and c+. In the particular case being considered, the inertia 
coefficient would be 0-95 for a free surface without gravity, and 1-05 for 
a rigid surface. Its variation with velocity can be seen from the curve A,, 
which gives the resistance R,. The coefficient p begins with the value 
0-95, rises to a maximum of about 1-07 near c/(gf )* equal to 0-4, falls to 
a minimum of 0-78 near c/(gf)* = 1-4, and then rises towards the value 
0-95 with increasing velocity. 

Similar calculations were made for the case y/g = 1-276, for which the 
wave resistance is shown in the curve B,. The curve for R, in this case 
is not shown in the diagram, because on the scale its magnitude would be 
nine times that of A,. However, the curve, in its relation to B,, is of the 
same type as in the case of the curves A, and A,, but with less variation 
in the inertia coefficient; this coefficient begins at 0-95, rises to a maximum 
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of 0-975 near c/(gf )! 1, falls to a minimum of 0-91 near ¢/(gf )* = 2-5, 
and then rises gradually towards 0-95. It is of special interest to notice 
that while the coefficient begins with the free surface value its rise towards 
the rigid surface value occurs before the wave resistance R, has become 
appreciable. A few calculations were also made for a very small accelera- 
tion, with y/g = 0-035, to confirm the general trend of the variation; in 
this case, the coefficient p has risen to a value of 1-05 at about c/(gf )* = 0-2. 

teferring to (37), some of the approximate values found in these 
calculations are given for reference. 

For y/9 0-1418, the quantity 32/87) has the value —0-52, —1-2, —1-1, 
0-92, 1-66, 0-4 for 6B equal to 40, 10, 4, 1, 0-5, 0-25 respectively. For 
yg = 1-276, the values of 3248) are —0-04, —0-24, 0-3, 0-4 for B equal to 
10, 1, 0-25, 0-125 respectively. 

The motion which has been examined in detail is uniform acceleration 
starting from rest. Similar calculations could be made for other cases of 
variable velocity, in particular for motion with uniform acceleration with 
a given initial velocity. In the latter case the results are not likely to be 
much different in general character; it appears that in any case the initial 
value of Rk, would be the inertia resistance for a free surface without 
gravity, and its subsequent variation would be similar to that shown by 


the present calculations. 
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THE ELECTROSTATIC FIELD OF TWO EQUAL 
CIRCULAR CO-AXIAL CONDUCTING DISKS 
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[Received 31 May 1949] 
SUMMARY 


In the earliest discussion of this problem Nicholson (1) expressed the potential as 
a series of spheroidal harmonics with coefficients satisfying an infinite system of 
linear equations, and gave a formula for an explicit solution; but this formula 
appears to be meaningless and its derivation to contain serious errors. 

In the present paper, starting tentatively from Nicholson’s infinite system of 
linear equations, a much simpler, though still implicit, specification of the potential 
is developed ; this involves a Fredholm integral equation the existence and unique. 
ness of whose solution are deducible from standard theory. The specification s 
obtained for the potential is shown rigorously to satisfy the differential equation and 
boundary conditions of the electrostatic problem. 

The Neumann series of the integral equation is shown to converge to its solution, 
so that the potential, and other aspects of the field, can be explicitly formulated 
and thus computed. 

The errors in Nicholson’s process of solving his system of equations are exhibited 
in detail, and it is concluded that attempts to carry through that process without 
error cannot lead to an explicit solution. 


THE subject of this paper has been discussed by Nicholson (1) and by 
Nomura (2); but on account of certain defects mentioned below and 
discussed subsequently, there seems to be room for further treatment of 
it. First of all, however, I wish to acknowledge my great indebtedness to 
Nicholson’s remarkable work; and to state that, although I criticize some 
parts of it, other parts are of inestimable value. 

Nicholson’s chief aim is to obtain exact formulae for the potential at 
any point. But his main result (N 91)+ contains an integral which is not 
convergent;{ and no indication is given, or has been found, as to how this 
is to be understood. Nomura aims chiefly at numerical calculation of 
capacities and total forces; but, although his results appear plausible, the 
convergence of many of his complicated expressions is not obvious, and 
is nowhere examined. 

In this paper I establish a new expression for the potential, involving 
the solution of an integral equation of well-known type, much simpler 
than that considered by Nicholson. This expression for the potential is 

+ This refers to equation (91) of Nicholson’s paper. Most of the equations so quoted are 


reprinted later in the present paper; see pp. 431-3 and 444-8. 
t See p. 444, below. 


(Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 4 (1949)] 
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shown to have an unambiguous meaning, and is verified to satisfy the 
differential equation and boundary conditions of the electrical problem. 
Further, it is accessible to computation because the Neumann series 
iteration) solution of the integral equation converges in all cases, and, 
indeed, at a known minimum rate. 

A detailed critical study of some extracts from Nicholson’s paper seems 
to be necessary in order to show the need for, and the advantages of, the 
present re-formulation of the problem. In this part of the paper I show 
that his results are indeed anomalous; and seek to explain why they are, 
ind to point out certain other errors of principle as well. I also show how 
, constructive attempt to carry through his methods of solution, without 
their errors, leads unfortunately to self-contradiction. It is almost im- 
possible to be brief in these criticisms, as the controversial situations 
arise late in Nicholson’s work and cannot easily be detached from the main 
body of it. 

As regards his infinite system of linear equations (N 19), Nicholson’s 
york probably admits a rigorous justification, even though he uses 
transcendental differential operations on infinite series (N 3) to establish 
his transformation formula (N 6), which is the corner-stone of this part 

fhis work. Nomura also obtains this same system of equations (see his 
7), (10), (12), (25), and (N 7)), but independently of Nicholson’s work, 
vhich he considers (3) to be too complicated to serve as a basis for com- 
putation. Nomura proceeds at once to numerical solution of the infinite 
system of linear equations by iteration, implicitly trusting that the process 
s convergent; this is the basis of all his numerical results. Nicholson, 
however, retains hope of an exact solution, and transforms the equations 
nto a variety of forms, (N 29), (N 43), (N 54), (N 61), leading finally to his 
momalous solution (N 91). It is these areas of Nicholson’s work which 
| have to criticize, in spite of their undoubted interest and partial validity. 

My interest in this work arose in the first place in consequence of a 
juery by Dr. V. D. Hopper, Dept. of Physics, University of Melbourne, 


in connexion with his measurements of the charge of an electron (4). 


Guide to the new developments 

Two leading cases of the problem are considered. They are: to specify 
the field due to two equal circular co-axial conducting disks at equal, and at 
equal and opposite, potentials, the potential at infinity being taken as zero. 
Nicholson and Nomura give most of their attention to these cases, since the 
lield due to the same disks at any given potentials is deducible from them 
ierely by superposition. They may clearly be called the cases of equally 
harged disks and of oppositely charged disks. 
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I begin by tentatively accepting Nicholson’s work as far as his specifica. 
tion of the potential by a series (N 20) in which the coefficients satisfy an 
infinite system of linear equations (N 19). I convert this to an integral 
formula (1) for the potential, involving a function known only to satisfy 
a certain integral equation (2). From here on I disregard the origins of 
this new specification, and endeavour to prove independently that it 
determines uniquely a potential function satisfying all the requirements 
of the problem. This analysis appears to me to be a more adequate 
justification of the result than checking the validity of the steps which 
led to it. 

The potential is expressible explicitly only if the solution of the integral 
equation is so expressible. It happens that the Neumann series gives such 
an expression of the solution, for it is convergent whatever the relative 
spacing of the disks. For computation one may start with any function 
thought to approximate to the solution; then the closeness of approxima- 
tion of its successive iterates is readily predetermined. 

The results established are as follows, the upper sign referring to th 
case of equally charged disks and the lower to that of oppositely charged 
disks. 
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Fic. 1. Dimensionless coordinates p, ¢, ¢’. Spacing parameter x. 
Upper disk at potential }. Lower disk at potential +). 


THEOREM 1. In the two leading cases described above, the potential at any 
point (p,¢,¢'), specified as in Figs. 1 and 2 by its distance r = pa from the 
axis of the disks and its axial distances z = Ca and z' = C'a from their 


planes, is 
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ifica here each square root has positive real part, and f(t) is the solution of the 
ry al itegral equation 
egy 1 
Le k 
tis f(a) . f(t) dt i (il c2#< th (2) 
s | «< d )2 
ns A 
ab i here « is the spa ing parameter | Fig. l). 
le! . . 
THEOREM 2. For every positive x, equation (2) has a continuous solution, 
jua . . - = eee i 
1 nd no other solution: it 1s real and even, and is specifiable by the Neumann 
hi : . 
res (Séé also Le mma 5S) 
1 
€01 f(x) 1+ ¥ (+1)" | K,(x,t) dt, (3) 
SUC ’ 1 Zs 


ATL | hore the iterated kernels K (v.t) are give n by 


CT | 1 
; : . l K : ~~ fe, 
ima- J K, (2x, t) ™ K,, (x,t) = | K,, -1(#, 8) K,(s, t) ds. 
7 K* ~—t)* 
~% 
0 tI 7 - sai ‘ 
[HEOREM 3. The capacity of each disk in the two cases is 
- y O 


1 


+ 


f(t) dt, (4) 


a 


. 


l 
nd the components of the field at all points not on the disks are given by the 


yproprvate jor mal differentiations of (1). 


Origins of the new formulae 

We now survey the temporary foundation, already constructed by 
Nicholson, on which the new formulae are to be built; see his paper, 
pp. 318-20 (and also 308-18). After some minor adjustments to re-orient 
t suitably, it is set out in the following paragraph. 

When the upper and lower disks are at potentials V, and +K, the 
potential V at any point whose spheroidal coordinates are (u,7) with 
respect to the upper disk and (’, 7’) with respect to the lower (see Fig. 2 


ind below) is, in terms of Legendre functions, 


, Moro ia ‘ “a? ee 
J a "beads | (1)V5,(¢n)A Py, (u Vo, (ty hs (N 17, N 20) 
where the coefficients a,,, are real and satisfy the following infinite system 
of linear equations with real coefficients: 
Ayn t(4m+1) > Kea, 6. (se = 0,1, 2....). (N 19) 


n m 
n=0 


Here 5, = 1, 5, 0 if m #0; and K™ (= K®(c/a)) are the coefficients 
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in the following transformation formula from spheroidal coordinates (y’, »' 


to spheroidal coordinates (yu, 7): 


imP (w’)Qn,(in’) = 7 > i"(n+-4)K2, P,(u)P, (in) (Ne 


m 
n=0 
for points (u, 7) suitably near the upper disk. 

Some explanations are desirable, particularly regarding the systems of 
coordinates. The notation differs from Nicholson’s only in using it 
place of his ¢, and the more usual Legendre function of the second kind 
@,,(8) instead of . . 
oe 9,(8) = i"1Q,,(i8). (5 

The upper disk, specified in cylindrical polar coordinates (r,6,z) by 


r <a and z = 0, is taken as ‘focal disk’ 7 = 0 of spheroidal coordinates 
(u, n); these are such that —1 <p <1, 7 > 0, and 











r = ay{(1l—p?)(1+ 7?)}, 2 = apy. (6 
ei ESSE 2x -const. pr const 
oe ~~ _ 1 (pos) 
\ ae een » 7 
Oe a: 
“et A? —_ =, "a 
x en ~+# Xx \ 
Hte* lL  z=0 7-0 } 1 7X\ \ 
ws | - 1 \ Wii r 
“= Pt ee) 
a ; ae 
aN fe ae lar ee ae P. 4 ?~ 
~ 
.™* Sigg ng O . 
/ Thi = agit c 
Re en -const 
zh (neg) 
2%=0 ” 
o' ~ oe 
<——-¢ >i< 2 >| 





Fic. 2. Nicholson’s coordinates. Upper disk at potential %. 
Lower disk at potential +-V. 


When the disks are charged, the potential due only to the distribution 
of charge on the upper disk is assumed to be the (formally) harmoni 


function OV. 
“v0 


7 


Ay Fy (H)In(7). (7 
n=0 
The coefficients a, are supposed real since q,,(7) is real for real ». Further, 
since replacement of 1» by —y merely replaces z by —z in (6), and so leaves 
the potential unaltered, a, is supposed to be zero for odd n. 
A similar system of coordinates (7’,z’) and (y’, 7’) is associated with the 
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be, lower disk. Since the disks are identically or oppositely charged, the 
potential at (’, »’) due to the lower disk only is, by (7), 


(N 4 2, , . : 
N 10'S a, P.(u")dn(n')- (8) 


Adding (7) and (8) we obtain (N 20). 
The boundary condition at infinity, 

V = Of(r?+2?)-4} = O(1/n) as > +00, (9) 
is already formally satisfied by (N 20). The boundary condition on the 
upper disk, which by the symmetry implies that on the lower disk also, is 

V>VW as n>-+0. (10) 
inate \pplying (N 6) to express V in terms of p and y for sufficiently small 





| positive 7, and formally equating coefficients of P,,(), Nicholson expresses 
. this boundary condition as the system of equations (N 19). 
| Most of the rest of Nicholson’s paper is concerned with solving equations 


N 19); for the present we do not consider this. 


Nicholson uses the following explicit formula for the A7,: 
K” = Ke) = [ ep sle yale), (N 7) 
\ wer 
0 
which shows that K”, = K For the new developments the following 
expression is more appropriate 
oem 
K" (x) = - Q,, (ix —t)P, (t) dt. (N 15) 
7 J 
\ proof of the transformation formula (N 6) with this form of K®, is needed 
later and is given in Lemma 2, partly because of the uncertainty caused 


by misprints, transcendental differential operations, and lack of rigour 
n Nicholson’s paper. The simple condition of validity found in Lemma 2 
does not appear in Nicholson’s account; this is not surprising in view of his 


methods at that stage. 


Development of the new formulae 


ut On the basis set forth at the beginning of the preceding section we now 
nol obtain the formulae of Theorem 1. In the following section they will be 
justified independently and rigorously. 
LEMMA l. Atall points not on the upper disk, that is, for which 
l LU 1 and yn > Y, 
ther 1 
AA TES ‘ a F P. (¢ dt 
si Py, ()Q> v1) : =" ) 5 ’ (11) 
; ’ 2 (2+ (2+ iat)?) 
vd (2-2 ; 
' —1 
h the 


there r and z are given by (6) and the square root has always positive real part. 
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First suppose z > 0, so that 0< p< 1 and »>0. By an addition 
theorem in Legendre functions (5), 


7 


. ] Z—r cos 
Py (4) Qon(t7) a. 9. Qxn( = — *« dd (12) 
] aP,, (t) 


dt 


I 





as { 


—-r -1 


4a | iz—r cos d—at 


by Neumann’s integral (6) for Q,,,( The onder of integration can be 
inverted, for the sen is ews since z ~ 0; thus 


wal 


7 


F,. ¢ 9 (2 = | PF... t) dt — 
On| L)Qsn( n) aai (z-+iat 


dd 
)+ircos¢ 


a [ P,,, (t) dt (13 
V 


{(z-+-iat)?— (ir)? 


by Jacobi’s integral (7). The square root is a continuous function of f, 
and takes the positive real value when ¢ = 0; for the inner integral which 
gives rise to it has both these properties. By comparing the arguments 
of (z+-2at)?+-r? and (z+-iat)? it is easily deduced that the square root has 
positive real part. 

Now suppose z < 0, so that 0 > » > —land 7» > 0. Then 


on(H)@on(?n) = (—1)?"Pan(—#)Qon(tm) 


Q.,{ =e" *) dé 


a 





by (12), replacing » by —y in (6). Then, by the same argument as before, 
(13) follows, with z replaced by —z which is positive. 

The integral on the right of (13) is unaltered if z is replaced by —z and 
t by —t, because P,,,(¢) is even; thus (11) is obtained for z < 0 also. 

It remains only to establish (11) for » = 0 and y > 0, that is, for points 
in the plane z = 0 of the disk but not on it. This ellenwa from (11) with 


p #0 by continuity. For the left member is a continuous function of 


(uw, ») at such a point, and hence of z; while the integral on the right is 
also a continuous function of z because its integrand is a continuous 
function of (z,t), for r > a, since 


9 


Ri{r?+- (z+-iat)*} = r?+-22—a2#? > r?—a? > 0. 
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LemMa 2. If x = c/a and kK" = K"(k) is defined by 


1 
gmintl ; — 
K” («) Q,,(ix—t)P, (t) dt, (N 15) 


= 


then the transformation formula (N 6) holds at all points (see Fig. 2 and (6)) 


it wha h 0 Y) t. 
Since x > 0, Q,,(ix—w) is a regular function of w = u+iv in the half- 
plane v x; it has therefore (8) the expansion 
4 
Q,,(ix—w) = > k” P(w), 
n=0 
miformly convergent in any closed ellipse with foci +1 wholly inside 
x, and, since (N 15) is to define Sl 
1 
: ; (n-+43)7 ,- 
oe (n+4) | Q,,(«—u)P,(u) du = ~—=— K®., (14) 
2 2 i gmt+n+l 
= 
By the addition theorem (12) already quoted, 
; 7 lc iz’—r’ cosd _ 
P(u OP m (en ) = Cn - dd (15) 
ertainly if 0 < p’ < land 2 0, (r’,z’) being related to (u’, 7’) as in (6). 
Writing 2’ ( Ae r, and c/a XK, 
( 12" r’cosd ° r cos d- iz 
Lm Qn —— 
a ] a 
S jn p re) 
~~ 2 ? 
a 
n 0 
ind this series is uniformly convergent in —z < ¢ <7 provided the 
ellipse with foci +1 through the points (+r—iz)/a lies wholly inside 


k; this requirement becomes, after a little calculation, just » < x. 
Then, by (15), 


; kn r cos d—iz 
Pr(t' Qmn(in’) = > 5" | P, dp 


a 


“a7 
n=0 oe 


> h(—"P, (w)P,, (in), 


using the well-known addition theorem in Legendre polynomials. Now 
N 6) follows using (14). 
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Using Lemma 1 we now obtain (1) of Theorem 1 from (N 


as follows: 


: 1 1 
pat Sonim f Pate f Pinar 
vir? 


a i Ly f (z-+-iat)?} ra v{r’2+-(2’4 iat)?} | 
hf 
; =, 
= - [ ¥ 2 “a “4\2) rm (2 vr . | - As, "Py, (t) dt, 
7 Y \/{p2 + (C-+12t)*5 Nie" +6 4 it)?\J bray 
-1 
using the notation of Theorem 1. We thus have (1) if we define 
f(t) > Aen i>" P,,(t), (16 
n 0 


necessarily an even real function. We need not consider in what sense the 
series (16) is to be understood if it is not convergent ; for the results of the 


present formal work are rigorously verified later. 


Using Lemma 2 we also obtain (2) of Theorem 1 from (N 19) formally, 


as follows: 


va 


— y2m 2n 
Sn, — Bay, +(4m i 1) S Asn ‘| Qom (tk t)P,,,(t) dt 
n=0 
9 1 1 < 
y2m - 
— (4m ; | Vom (UK t) Ss As, P, ( ) dt 
5s . &, 
1 
y2m ae 
(4m+1)— ] Qoplix—Of(O) dt 
oe 


l —f(x) D2 (8,,.— Bom)t?™ Pom( v) 
m=0 
i< 
eae / (4m r om\ | Qom(tK t)f( t) dt 
cis L—" 


3 


a f(t nay (4m-+1)P,,,(x)Qo,, (ik —). 


m n=0 
oj 


Now (9) for —l1 <a <1, 


> (2n+1)Q, (ix—t)P, (x) = ae 
n=0 uK—t—Z2x 


20) formally, 
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ally, # since the ellipse with foci +-1 through ix«—# contains —1 < x < 1 in its 











interior. Also, replacing x by x, 
; ] 
| ¥ (—1)"(2n+1)Q,(ix—t)P, (2) = ———, 
| n=0 u—t+-x 
so that. for l a l. 
] ] ] 
S (4m L)Q (2K i Fe, (x) = ——— 
0 2\on—t+a °° u—t—a 
Thus, fo1 l a a 
1 
i ] ] 
l (x) — -) f(t) dt 
20 | u—ttx u«—t JS 
l 
l 
F i f l ] ‘ 
= ————) 
if +] 2a | ie t+a2 «u«-+t ys 
l 


replacing ¢ by —¢ in the second term since f(t) is even. This immediately 


Proof of Theorem 1 

It might be easy to justify the three formal integrations of series 
involved in the preceding section, under a suitable assumption as to the 
order of magnitude of the unknown a,,; we should then have a valid 
proof that (N 19) and (N 20) imply (1) and (2). But (N19) and (N 20) 
were found by Nicholson using a number of formal steps needing justifica- 
tion (mainly double-limit inversions associated with (7)—(10)), and in view 
ff the defective result given later in his paper it seems preferable to 
dispense with them except as scaffolding leading to the formulae of 
Theorem 1. 

We accordingly establish Theorem 1 independently, by verifying that 
Laplace’s equation and the boundary conditions of the electrical problem 


are satisfied by (1) if f(¢) is a continuous, real, and even function satisfying 
2). It is then necessary to prove that such a solution of (2) exists, and 
that it is also unique; this is done in the following section, where Theorem 2 


18 proved. 


Lemma 3. Jf f(t) is continuous, real, and even, and the real part of the 
square root is positive. then the real function of (r, 0,2) 


a [ f(t) dt 
T yr? (2-4 iat)*' 


« 
1 


S hay monic ¢ verywhe re exce pt on the disk Z V0, ys @. and is continuous. 


jor normal approach from either side. on this disk. 
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The Laplacian of this function can be expressed by means of the 
appropriate differentiations within the integral sign, provided the inte. 
grands so obtained are continuous functions of (7,z,¢). This is so provided 
no denominators vanish, that is, 


and this condition is satisfied everywhere except on the disk z = 0, r <4, 
Hence the stated harmonic property, since the Laplacian of the integrand 
is zero for each value of f. 

For the continuity on the disk, suppose first that 0 << r <a. Then 


r2+-(z+-2at)? J {(r?+-22—a#?)?+ 427070} 


~ 9 


z2tr2—q%y? if att? < 7, 
> r2—q?f?; 

|r2+-(z+-iat)?| = ./{(r?—z?—a#?)?+4 427} 

> 22+a%—r? if a*t? > r?, 
. a*t2— 72; 


1, 


//\ 


whatever real value z may have. Thus, for —1 <¢ < 


fo es \f@)| 


~ 1 
\|r?—a?¢?| 





and so, by the theorem of dominated convergence, as z > +0, 


1 1 

a f(t) dt a fi) dt (17) 

a J Vfr2+(z+iat)> a2 J (r?—at?)° 
at | 


This establishes the stated continuity on the disk, except at r = 0. The 
square root on the right of (17) is, of course, to be understood as the limit 
of that on the left; it is thus the positive root if t is between +-r/a, but 





+i/(@?—r?) if ¢ = +r/a. 


_ 


Since f(t) is real and even, the right side of (17) thus reduces to 


a ’ f(t) dt ] r r 
_ — an -- Ss d ° 18 
z } J(rP—a)’ «| tz ‘) ? 


—rja 0 


In the case r = 0, until now excluded, the integral in question is 


1 
a [ f(t)dt 
Ca | |z|+-cat 


at 








sine 
mT 





of the 
€ int 
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since the square root is to have positive real part and interchange of ¢ and 
_t makes no difference; thus the integral is 


z|—vat 


5 { ( 2 7 iat) aM dt 
—1 


1 
3 | —_ 249 f(t) dt 


z|“+at 


“1 
> J(9) 
as > 0, by the usual argument with positive ‘ ,ular kernels like those 


of Fejér and Poisson. This shows that (18) still gives at r = 0 the correct 
interpretation of the left member of (17) to make it continuous on the disk, 


Lemma 4. If f(t) is continuous, real and even, and satisfies the integral 
equation (2), then (1) takes the value V, at every point of the upper disk and 
the value +V, at every point of the lower disk ; and at infinity (1) is O{(r?+-z?)-4}. 


Asz—>0 with 0 <r<a, 
¥ ] _ 1 
Hf_ fd  _ he [ f(t) dt 
aw) v{p?t+(C+it) oa J fr?+(2+1at)} 
1 _ 
ai A [ socos4) dd (19) 
0 


by Lemma 3, using (17) and (18) in particular. 


The other term in (1), omitting the prefix +, 


it f(t) dt Ya f(t) dt 


——- 
_ 


am ) V{p?+-(C’+2t)?} TT J {r?+ (2'+-2at)*} 


1 om 


is certainly continuous wherever 2’ + 0, since it is harmonic by Lemma 3. 


Hence its limit as z > 0 is given (see Fig. 2) by putting 
¢, (’ =2'/a = c/a = «. 


Using Jacobi’s integral as in (13), the square root having positive real 
part as assumed in Theorem 1, this term becomes 
l 7 4 
V . ” d 
: f(t) dt ————— « (20) 
r J k+it+ip cos ¢ 


l 0 
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We may replace either ¢ by —t or ¢ by 7—4; since f(t) is real and even. 
this gives the two equal conjugates 





1 7 
hf _. 
bi hed ee om 
—j} 0 
a r id 
_% It vw 
73 [fo J «-+(t—pcos¢)? 
—j 0 


their arithmetic mean, which is equal to (20). Also the order of integration 
may be inverted, since the integrand is continuous, giving for (20): 


1 
Re K 
Id | - . ' 
| Ft (peoag— ye) “ 


Thus, as z > 0 with 0 <r <a, the limit of (1) is 


a 


K 


K- (pcos dapat) so iat 





71 7 
2 


« 


0 ; —j 


4 | re, 
om | | flo coss) + 


Now this integrand is simply 1, by (2); for p = r/a < 1, and consequently 
—1 < peos¢d < 1. The whole expression thus reduces to J), as stated. 
Similarly on the lower disk (1) takes the value +), by its symmetrical 
relation with the two disks. 
For the behaviour at infinity it is enough to show that the integral of 
Lemma 3 has the stated order of magnitude. Since 


|r?+-(z+2at)?| > Rif{r?+ (z+ cat)?! > r?+-22—a?, 


the integral of Lemma 3 
] a ft 


= f(t)| dt, 


< V(r?-+-22—a?) 7 


. 


1 
which has the stated order of magnitude as ,/(r?-+-2?) > 0. 


Proof of Theorem 2 

In Lemmas 3 and 4 we have established that (1) satisfies the differential 
equation and boundary conditions of the electrical problem if f(t) is con- 
tinuous, real, and even, and satisfies the integral equation (2). 

We now reduce these requirements to the single one that f(¢) should 
satisfy (2), by showing that it is then necessarily continuous, real, and even. 
This is really needed to complete the proof of Theorem 1 in the precise 
form in which it is stated; and it also gives what is more conveniently 
stated explicitly in Theorem 2. 
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Lemma 5. If @(t) is integrable, x and d are positive, and p = «+A, then 


1 x 1 
- if ; . . 
| —_. p(t) dt — A ds = -D(t) dt. 
J pet (x—t)° A?+-(a—s)* K?+-(s—t)* 
-1 


1 
The order of integration on the right can be inverted, since this gives 
an absolutely convergent integral, namely 
1 
M(t) dt 


J ‘ 
The inner integrand has simple poles in the upper half-plane at s = « +-2A 
Hence evaluating the inner 


and 8s t--ix, since A and « are positive. 


integral, the right side can be written 


l l l ] | 
iN e—t+ia—ix  x—t+iA+ik)’ 


1 
| @(t) dt} 


21 lan—t—ice—iA 2—t—m- 


1 


which reduces to the left side as stated. 


Lemma 6. If @(t) is integrable and satisfies, for —1 <x <1, the homo- 


reneous inte gral equation 


P(x) -D(t) dt = 9, 


then (ft) is ide ntically zero. 
Suppose @(¢) is a complex integrable function satisfying this equation. 


Then 
K | (t) | dt. 


M(x) mY 


T K2-+ (x 


Substituting from this inequality into itself, and using Lemma 5, 


_ 


D(a . |D(t)| dt 





x?-+-(s—t) 
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By repetitions of this process, Ch 
1 1 , 
it 2M 1 once 
O(x)| < : , |@(t)| dt < —— | |@(t)| dt, mere 
am) (2"«)?+(¢—t)? 2"«a J 20 
ne = diffie 
which tends to zero as n 00; hence the result. doub 
To 
Lemma 7. There is a function f(t) satisfying 
1 
‘ a i K - 
(x) + - - f(t) dt = 1 D) 
f(x) Ee) bla pel (2) 
1 
for —1 <2 <1, and tt is unique, continuous, real, and even. 
Since the homogeneous equation has only the one solution by Lemma 6, 
it follows from the Fredholm alternative-theorem that there is a solution 
of (2), and that it is unique. 
] K 
Write = K(x, t), 23 
a K?+(x—t)? i, 0), - 
1 | F 
then f(x’)—f(x) | | {K(a’,t)— K(x, t)} f(t) dt} — 
= this 
= 1 Si 
bnd | K(x’, t)—K (a, t)| | | f(t)| dt; | ieee 
; Li 
this tends to zero as |x’—x| > 0, and so f(x) is continuous. 
The imaginary part of f(t) satisfies the homogeneous equation (22) and | 
is therefore zero by Lemma 6; thus f(é) is real. 
Replacing x by —wx and ¢ by —#, it appears that (2) is satisfied by f(—z) | the : 


also. Hence f(—x)—/(x) satisfies the homogeneous equation (22), and so, | to it 
by Lemma 6, f(x) is even. tern 


Lema 8. If f(a) 1s any continuous function and 


K 


Fr(@) — KL (g— petal) dt, ) A 





then the solution of (2) is lim f,,(x); in the notation of Theorem 2, 
no 


n—1 , 2 
j{«) = 14 > (1) | K,(a, t) dt+-(--1)" K,,(«, t) fo(t) it), Thr 
v 1 . . 
. —" infir 
— » DD vicianiianes 
and bnd | f(z) —f,(x)| < (: arctan } bnd | f(~)—fo(x)|. (24) 
7 K 


With fo(x) = 1, these formulae relate to the Newmann series (3) of Theorem 2 
and its remainder after n-+-1 terms. this 





a 6, 


ion 
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Clearly f,,(x) is uniquely determined by the stated recurrence relation, 
once f,(x) is assigned; so the stated expansion of f(a) can be established 
merely by verifying that it satisfies the recurrence relation. There is no 
difficulty in doing this; the order of integration can be inverted in all 
double integrals encountered since their integrands are continuous. 


To prove (24) we have, since f(x) is the solution of (2), 


Lf K " ; 
J(®%)—Jn(*) 3 4a pet Sn—i(t)} at 
1 
Z 
K 
f(x) —f..(x) bnd | f(é) (t dt 
' In—1\¢) mJ «?+(x—t)? 
1 


l{ l—zx . 1+2) 
bnd | f(¢)—f,,-,(¢)| — arctan + arctan — 


77 \ K K 
2 ] 
bnd | f()—f,,-(¢)| —arctan -, 
7 K 
since the penultimate line is greatest at x = 0. Repeated applications of 


this result give (24). 
Since arctan(1/«) < 4x, (24) shows that lim/,,(x) is f(a), the solution of 
2), as stated 
Let fp(x) 1 (—1 <x <1); then 
n ; 
Fin(®) 1+ > (1) | K,(x,t) dt, 


— 


y=] . 1 


the sum to n+ 1 terms of the Neumann series (3). Hence f(x) is the sum 
to infinity of (3), as stated in Theorem 2; and the remainder after n+-1 


2 1\" 
arctan . 
TT K 


As in the proof of Theorem 1, the potential due to the upper disk only is 


terms is of the order of 


Proof of Theorem 3 


V, l | fiat 


mJ y{p?+(6+it)?} 


Thus the total charge on the upper disk is the limit, as (r?+-2?)? tends to 
infinity ot 


1 

et (t) dt : es 21 9ale 

(2420 [ afd | f(t) dtl -0(" tiene )|; 
T r-+ (z-+-2at)“}? 7 7 \ r-+-z2° 

—1 —1 


this gives (4) as the capacity of the disk, since V, is its potential. 
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The components of the field due to both disks together, at all points not 
on either, are given by the appropriate formal differentiations of (1), as 
in Lemma 3. 


Criticism of Nicholson’s result 

The final conclusion of Nicholson’s treatment of the case of equally 
charged disks is a set of expressions (N 90-94) for the potential at} (r, z, 2’), 
which we can summarize as 


fe-iSite-IS18,J. ( pt)D(t) dt. (N 90-94) 





Here ®(x) embodies the coefficients a,,, of (N 20) in the form 


D(x) : S Asn Jon t (x )am - (25) 
0 


n 


and it is specified in a form equivalent to 


P(2) E ) as anh(sz Palas xs ds. (N 88) 


tanh(7/2«) 


The question of what méaning can be attached to these formulae arises, 
In the ordinary sense (N 88) is meaningless; for 





tanhas 1 log| anh As+tanhaA A+ O(e-2) (26) 
tanh a ~ ~|tanh As—tanh A} 
as s > +, so that (N 88) fails to converge at the upper terminal. 
In the sense of convergence in mean pth power, where p > 1, (N 88) is 


also meaningless; for 
x 


| sin as ds 
0 


is not convergent in mean, hence neither is (N 88) by (26), because 








0) 28 S p ox ' . 
we i ae - leos Sa—cos 2Sa|” 
| sin ws ds sinasds| dx | dx 
| a | 
0 0 0 0 
x 
- lcos t— cos 2¢|? 
Sp-1 | : dt 
t 





which tends to infinity as S > co if p > 1, and is already divergent if p = 1. 
In the sense of Abel summability (N 88) has a meaning; but the potential 


ft See Figs. 1 and 2 for the notation: r = pa, z Ca, z’ C’a. 








V th 
this 


For, 


Out 

Tl 
mail 
equi 
orig! 
step 
case 
case 


(a 


is il 
stag 


(N: 


an 


Thi 


equ 


obt 


to 


94) 


OF 


( 


Ss) 


(96) 


5) 1S 











TWO EQUAL CIRCULAR CO-AXIAL CONDUCTING DISKS 445 


V then given by (N 90-94) is meaningless at every point of space because 


this integral behaves like 


For, interpreting (N 88) by Abel summation, and using (26), 


? , ; tanh As 

(427)'D(x) lim | ( Qa 

" S5>+0. tanhaA 
Uv 


}sinars ds = =4 O(1). 
Outline of Nicholson’s procedure 

This outline is given in order to provide a basis for criticism. It refers 
mainly to Nicholson’s attempt to solve the infinite system of linear 
equations (N 19), and condenses drastically his exposition which, in the 
original, occupies pages 320-48; but it seeks to preserve the principal 
steps in sufficient detail for subsequent discussion. It refers only to the 
case of equally charged disks, since Nicholson discusses fully only that 
case. Many obvious misprints have been altered without comment. 


x 


(a) The function D(a > Gnd y(x)a-t (25) 


is introduced immediately, although often written as f(x)a-* in the early 

stages (e.g. (N 42)). The problem of determining the coefficients a,,, of 

N 20) to satisfy , 

Ayam t+ (4m-+1) > Ke" Asn Sn (N 19) 
n 0 


is then converted into that of determining ®(x) to satisfy 


. 1D ij 
O(y) H(y, x)e-**@D(x) dx R —. (N 29. N 42) 
. 7 Y 


in integral equation whose kernel is related to, but differs from, 


l(sin(a+y) | sin(a—y) s 
H(x,y) = -| #4. y)\. (N 28) 
w\| x+y x—y | 
This function is the kernel of the ‘half-range’ form of Bateman’s integral 
equation (10) 
, l f sin(y—ax i 
WY Y ) (x) dx (27) 
a y¥—2z 
obtained by the obvious transformation when 7(2x) is even. 
(6) The functions J,, v)v-* are all solutions of (27), so it is reasonable 


to assume that @(2) is also a solution; then by means of 


@(y) 2 sin y | Hiy.2)(e)— |- sin x| lie 
J \ 7 x j 


N77 Y 
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(N 42) is converted into an integral equation of the first kind, 


} iy 1 ¢ “)(2)— [= sin) (x,y) dx = 0. (N 43, N 50) 
2 ra 
(c) Applying, in effect, Parseval’s theorem for Fourier cosine trans- 
forms, (N 50) becomes 
1 


5 F 2 sinx 

| cos ya da ly + e-*2) D(x) — loos ve dx = 0. 
} J \ 7 «x | 
0 0 





Equating the inner integral to zero, this is certainly satisfied if 


x 


| (1+e-**)@(x)cos ax dx = Jz (O'<= a =< FE). (N 52) 


0 
Now it is convenient to write 2q in place of «, so that 
1+e-*t = 2e-%* cosh qx. 


(d) By formally inverting a transcendental differential operator 
cos(qe/éx), (N 52) is replaced by 


or 
ad 


= 8 
=) 9 ha . 
qc 2 { sina 
2 | e-@O(x)cos ax dx sec d - | cos ax dx 
OX 7 x 
0 


0 
x 
2 f sing cosax vn 
re | Y dx. (N 54) 
a) «x coshgqx 
0 


An alternative argument, not involving such an operator, is given subse- 
quently, in (N 62)-(N 63); it will be discussed presently. 

We now quote two formulae (11) repeatedly used by Nicholson: If n is 
a non-negative integer and x and y are real, y positive, then 


% 9\1 
| elle+iun .(t)t-# dt (=)Pimr0,(e-+ iy), (28) 
, 7m 
, 2\}. ; 
| ett], (t)t-* dt ( ins AQ) (a+ 02); (29) 


remembering that 
Q,(7+01) = Q,(x)—4mP(x~) (—l<a4< 1) (30) 


while outside —1 < x < 1 we must replace P,(x) by zero. 
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(ec) Substituting (25) in (N 54) and integrating term by term, the left 
side of (N 54) becomes 


~ ly, f i x21 p—(Q tia)e J. s(x )ar t dx 


> aon{-)° 2" Qon(tgt+-a)+Qon(ig—a)} 
0 
provided Im(iqg-+-«) > 0 so that (28) is applicable; thus (N 54) becomes 


' 
j 


> Ao», yan Qo, (0g +) + Vs (2g x) 


" sin x sinh(2q¢+-a(q- x) jac +-sinh(2q-+-2(a¢ x) ay (N 59) 
J 2 sinh 2qx 


0 
and this is certainly satisfied if, for suitable values of p, 


x 


” 1 f sine sinh(2¢ x 
> do, 07"Qo, (tp) (24—p) dx. 


n=0 o. 
0 


(N 60) 


x sinh 2gx 


Putting ip = s, and evaluating the integral, this is rewritten 


tanh(sz/4q) — 
SY (—1)"a,, Q,,(8 ¢ ‘ N 61 
n=0 "Gan Sank) ec 4q) ati 


(f) Equation (N 61) is now regarded as analogous to a Fourier expan- 
sion of the function on the right in a series of orthogonal functions on 
-1<s < 1, and the coefficients a,,, are written down explicitly: 


. a 


' tanh(sz/4 a 


« « 


1 -1 
This completes Nicholson’s solution of equations (N 19). 

(g) The potential at any point is now known from (N 20) and (N 86); 
but it is also expressible in the more compact form (N 90-94) (already 
quoted), which can be made explicit without the use of (N 86). This is 
done by solving for ®(¢) by means of (N 61) as follows. Replacing n by 2n, 
the imaginary part of (29) gives, for any real s, 


> (1) "Gon Qon(8) = V(37) > Ae, | sinst J; 


2n 


0 
Using (N 61) on the left and integrating formally on the right, 


x 


/9 


2 é (Sar/ 4 eae i ‘ 
Qa nie... J | (t)sin st dt. (31) 
tanh(7/4q) N 


e 


0 
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By the inversion formula for Fourier integrals, then, 
x 


2 (2. /tanh(sz '4q)\.. ‘ 
~— Ris } it (Sante wt iad das 


allt 
0 
Criticism of Nicholson’s procedure 

The most natural check on the validity of the work is a direct verifica- 
tion that the differential equation and boundary conditions are satisfied 
by the final expressions (N 90-94) given for the potential in the case of 
equally charged disks. Lacking an interpretation of these expressions, 
this is out of the question. 

A less comprehensive check might be attempted, for instance to verify 
that equations (N 19) are satisfied by (N 86). This would be far from 
simple; it seems more worth while to study Nicholson’s process itself. 

Nicholson seeks to justify the step from (N 52) to (N 54) by an alternativ 
argument on p. 340, which we now reconstruct critically. 

First consider the equation 

: x 


F 2 f sine cosh px 
2 | e-% cosh pa O(x)cos ax dx J P. pe dt 





—— COS ax — 
7 x cosh qx 


0 r ; 
(N 62) 
supposing «a real and 0 < p < q to ensure convergence of both integrals. 
As p > q, (N 62) formally approaches (N 52), if 0 < «a < 1. In both sides 
of (N 62) substituting 
2 cosh px cos ax = cosh(p+-ia)a+cosh(p—ia)a, 
it appears that the two equations 
[ e~2 cosh(p+tia)a O(a) dx — 


* \V 
0 0 


— dx 


x 
] [ sin x cosh(p-+-ia)a 
ii cosh gx 


Da» 
together imply (N 62); but the converse is false. Thus (N 62) is implied, 
not reversibly, by 


i l sin x cosh Ba 
| e-2 cosh Ba D(x) dx = [ F ( 
e a) ( 7) J 
0 0 


meas . lx (32) 
2 x coshqx 
with B = ptia, 0 < p < qand a real. 

Now let pq; it is expected that (N 52) is implied, not necessarily 
reversibly, by (32) with 8 = g+ia and 0 <a <1. 

But, in effect, Nicholson puts B = ix in (32), and regards this as estab- 
lishing (N 54) (see his comment on (N 63) and (N64)). Actually he first 
applies to (32) the same treatment (e) as he has used on (N 54), and then 
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identifies the results of this treatment, (N 63) and (N 59), by putting 
ix; but the identification can clearly be driven back to (N 54) as 
ybove. 

Assuming the validity of the limit process from (N 62) to (N 52), the 
woper inference from the above argument is that (N 52) is implied, not 
necessarily re versibly, by (N 54) provided the latter holds for all values of x 
n the line-segment joining iq+-1. If O(x), independent of «, can be found 
to satisfy (N 54) for this range of values of «, then it necessarily satisfies 
N 52). 

We now reconstruct Nicholson’s subsequent steps with a view to salvaging 
the work as far as possible. 

Assuming the validity of the term-by-term integration involved, (N 54) 
s implied by (N 60), again not reversibly, provided the latter holds for 
0 = ig+a; that is, for all values of ip on the line joining +1 and on the 
line joining 2q7-+-1, to give « the values required above for (N 54). 

For the real values of ip thus required, (29) must be used in place of (28); 
so that, in (N 59) and (N 60), Q,,,(s) must be understood as Q,,,(s+02) if 
-l<s | (contrary to convention). For the real values of ip the term 
by term integration is more hazardous, because each term on the left of 
N 60) arises from a non-absolutely convergent integral like that in (29). 

Equations (N 60) and (N 61) are equivalent for the requisite values of 
ip, provided the Q,-function on the right of (N 61) is also understood in 
the above sense when —] s 1. A full verification of this is lengthy; 
but the imaginary part is easily checked, and this is what is needed for 
35) below 

Reviewing in summary the relations between (N 50), (N 52), (N 54), 
und (N 61), we see that (V 50) is satisfied by O(a) of the form (25) if (N 61) 
holds for all values of s on the line joining +1 and also on the line joining 
2gi+1; provided that Q,,,(s) is replaced by Q»,(8s+ Oi) on the former line, 


ind two double limit inversions involving the unknowns are valid. 
Thus the coefficients a,, are to be real numbers, independent of s, 
satisfying (N 61) in the following two forms simultaneously, for s between 
l and 1: 


S (—1)"ay, Qo, (8-01) = Qe tanh(sz 49). oi), (33) 
an ‘ 5 tanh(z 4q) 


(34) 


> ( 1) ‘Ay dV s 2q1) ot ee i). 


tanh(7z/4q) 


The imaginary part of (33) gives, by (30), 


1)", Pyn(8) 
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which requires (12) that a, = 1 and a,, = 0 (n 4 0). Substituting these 
values in the real part of (33), 


tanh(sz/4q) 
Qo(s) 0," t anh(a a (—i <4 < 2), (36) 


which is self-contradictory. Similarly (34) becomes self-contradictory, in 
both its real and its imaginary parts. 

Thus no values of the coefficients a,,, satisfy all requirements, because (33) 
and (34) lead to contradictions. 

Nicholson does not encounter this final impasse. This is because his 
treatment of (N 54) leads him to require that (N 61) should hold only for 
values of s on the line joining +1, not also on that joining 2gi+1; and to 
consider in effect only the real part of that equation. Thus (34) does not 
enter at all, nor does the contradiction arising out of (33). 

Nicholson goes on to use (N 61) in two ways, (f) and (g); in each case 
treating the Q,,,(s) as real functions. 

In (f), (N61) is mistaken for an expansion in orthogonal functions 
Q;,(s) on —1 <s <1, and formulae (N 86) for the coefficients a,, so 
deduced, They are, of course, erroneous; but that for a, is used in (N 89) 
and (N 112) as one of the conclusions of the paper; the same type of formula 
is also given in connexion with oppositely charged disks, in (N 132), 
(N 134), and (N 135). 

In (g), (N 61) is used to obtain the Fourier sine transform of ®(z), in 
order to deduce an expression, (N 88), for ®(2) itself. This too is erroneous, 
at any rate in the form presented, because the transform is known only 
for —1 < s < 1, since (31), like (N 61), is restricted to this range; this is 
of course insufficient to determine the function D(x). Nicholson assumes, 
without any apparent justification, that (31) still holds, in effect, outside 

l<s < 1; this cannot be, because the Fourier transform of an integrable 
function tends to zero as s > x, by the Riemann—Lebesgue theorem. It is 
because the left side of (31) tends to a non-zero limit as s > «, as in (26), 
that the integral (N 88) is divergent. 

One further point of criticism is that from (N 43) to (N 52) the integral 
equations used do not specify D(x) uniquely. This is particularly apparent 
in (N 52), which gives the Fourier cosine transform of (1-+-e-**)®(x) only 
for values of the argument less than 1; for other values it is left unspecified. 

This situation originates in the use of Bateman’s integral equation (27). 
Equation (N 43) is of the form 


f p(x y) dx = 0, (N 43) 


0 
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ind. defining us(a) to be even, can be written, like (27), 


] sin(y *) hi) dx = 0. 


1T Yy x 


« 


Nicholson, in § 14, mentions various particular solutions; the general 
solution in L? comprises a// Fourier transforms of even functions of L? 
vanishing in (—1, 1). 

A selection from this embarrassing variety of solutions is implicitly 
made after (N 52) by remembering that (2) is of the form (25). A further 
selection is also made on the criterion that as one disk is moved away the 
potential ne 


r the other should tend to the known potential due to a 
single disk. Analytically this criterion is expressed in the two forms: 


2 siny . 


D(y) J,(y)y-* as «>, 


NT Y 
a,—>1 and Ao, > 9 (n + 0) as x0, 
which appear in connexion with (N 29), (N 53), (N 58), (N 60), and (N 61). 
This procedure could be relied on if it were known that only one of the 
10st of solutions of (N 43) satisfies both the criteria of selection. 


REFERENCES 


1. J. W. Nicuotson, Phil. Trans. A 224 (1924), 303-69. 

2. Y. Nomura, Proc. Phys.-Math. Soc. Japan, 3rd ser. 23 (1941), 168-80. 

3. loc. cit., p. 168. 

4, V. D. Hopper, ‘The Electronic Charge and the Oil Drop Method’, Australian 
Journal of Scientific Research, ser. A, 1 (1948), 377. 

5. E. W. Hosson, Spherical and Ellipsoidal Harmonics (1931), 378. 

6. loc. cit., p. 63. 

7 loe. cit., p. 361. 

8. E. T. Copson, Theory of Functions (1935), 296-7. 

9, loc. cit., p. 296. 


10. G.H. Harpy, Proc. Lond. Math. Soc. (2)7 (1909), 445-72 ; and E. C. TrrcHMARsH, 
Fourier Integrals (1937), 349-51. 

11. G. N. Watson, Bessel Functions (1922), 387. 

12. A. Zyamunpb, Studia Math. 2 (1930), 143. 











TRANSFORMATION OF CERTAIN SERIES OCCURRING 
IN AERODYNAMIC INTERFERENCE CALCULATIONS 


By F. W. J. OLVER 
(The National Physical Laboratory, Teddington, Middlesex) 


[Received 18 January 1949] 


SUMMARY 
A transformation is established which converts a certain slowly convergent series 
into a rapidly convergent and easily computable form. Numerical values of the 
coefficients of the new series are given for some cases of practical importance. Thi 
method used may be extended in its application ; general forms of slowly convergent 
series which can be similarly transformed are stated. 


1. Introduction; 

In aerodynamic calculations of interference on lifting surfaces in rectangu- 
lar wind-tunnels (1, 2), the chordwise variation of induced upwash involves 
the computation of the following double series: 


x 


y > (—)"f,.(a—m)—f,(B -m)|, (1) 


| tom a 
n=1m=-—0 


F(a, B) 


/ 


y(y*+-2p?n?) 


where SY) ren 
n“(y>—- pn)? 
x, 8 are real parameters and yp is a positive parameter. In aerodynamic 


work y is the ratio of tunnel height to tunnel breadth. 

The series (1) is very slowly convergent and laborious to compute by 
direct summation. Brown (1) gives an estimate for the remainder after 
a finite number of terms, and though the use of this estimate considerably 
shortens the work, it remains laborious. 

In this paper the following transformation is established, 


F(a,B) = x(B)—x(a), (2) 
where / 

x(x) tr’xt > k,sin 2zsa, (3) 
s=1 

=. ‘ =. K, (27s 
re. 877s >» (—)"-1K,(27psn)+- 4 z (—)"-} i a (4) 
n=1 n 
n=1 


and K,(z) denotes the modified Bessel function (Watson, 3). 
The series (3), (4) are each rapidly convergent unless jz is very small. 
+ The author desires to acknowledge the collaboration of Mr. H. C. Garner who supplied 
the aerodynamic information in this section. 


[Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 4 (1949)] 
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[It will be seen that a transformation of this type can be established for 


any convergent series (1) with f,(y) of the more general form 


(r+) 
> 


Ws. (y, w)(y? +p?) 


where us,,(y, #) is a polynomial in y of fixed degree, and r is a fixed positive 
integer. The transformed series will always converge more rapidly than 


the original. Convergent series of the form 


tL ub(a—m, 2) = b(B—m, pu | 
} 


— {(a—m)?-+ pe? hrs {(B —m)?* tp? ir 


where %(y,) is a polynomial in y, may also be transformed into more 


rapidly convergent series by the same method. 


2. Transformation of the series 
Where fractional powers of complex quantities occur in the following 
analysis their principal values are to be understood. Write 


S) (x, B) > [f.(a—m)—f,(B—m)| (n > 1), 


so that F(x, B) > ( YS) (x, B). (5) 


Then by Cauchy’s residue theorem, 


S, (a, B) [_ $2) dz, (6) 
where 
[(a—z){(a—z)?+2u2n?! (B—z){(B—z)?+ 2u2n?}] cot mz 
= f(a—z)2+pent = {(B—z)®+ pn? Sar 
and C is the contour comprising the two straight lines |Imz C, 
0<e un (Fig. 1). The singularities of the function ¢(z) consist of 
simple poles at the points z= 0, +1, +2...., and branch points} at 


2=atipn, B+ipn. The cuts corresponding to these branch points are 
parallel to the imaginary axis (Fig. 1). 

In the domain |Imz c, (2) O(\z|-*) as z > 0, except between the 
pair of cuts in the upper half-plane, and also between the pair of cuts in 
the lower half-plane. In the former of these excepted regions ¢(z) > —n-, 
as z > 00; and in the latter 4(z) > n-*, as zoo. It follows that in (6) the 
contour C may be deformed, and 


S) (x, B) a x +-lim {9 (z) dz. ( 


1” pan 


~I 
~— 


7 It is assumed in Fig. 1 and in the analysis that « > B. That the final results remain 
true if a B follows by interchange of xX, B. 
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I’ comprises the aggregate of the four closed loop contours Tj, Ij, f,, I, Add 
(Fig. 1). the 
Consider the integrals round I, [,. Make the transformation § (0 
t = 1(z—a)/pn, 
/MAG. CUT CUT 
AX/S 
a+ip I 
c¢ 
can 
- La risn ly 
™ 
tc C, 
7 
REAL _AX/S 
-ic & He 
7 . 
of 
In 
CUT CUT 
Fic. 1. 
so that a—z = ipnt. Then by use of Cauchy’s theorem H 


(-1+) (1+) 
[ ¢(z) dz = — e [ - [ ti *) coth a7(unt+-ia) dt, (8) | A 
JN+h 2n} . J ja—#) 
. =@ qd v 
where g = p/un. 
The paths of integration on the right-hand side of (8) are illustrated in 
Fig. 2. Similarly 


i d(z) dz = ef ‘ +f ~ |) oth x( (unt+-7B) dt. (9) 
B+l =n 





| | | 
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Adding (8) and (9), and letting p >, equation (7) can be expressed in 
the form 


J (a, B) — (a 3) 4 “ sin ar(« B) x 
n* an 
(—1 (1 ) 
[ [ t(2 -t*) dt (10) 
: J |(l—@)! sinh a(unt-+-ia)sinh m(unt+ 7p)’ 


It is easily shown that 
cosech z(unt+-ia)cosech 7(unt-+-78) +-cosech m(unt—ia)cosech m(unt—ip) 


can be expanded in the form 


x 
4 cosec 7(a—f) > A,e-?7#8™, where A, = sin 27sa—sin 2zsf. 


1 
CUT a ; Ce wy 
’ 7 i] 


Fig. 2. 








Hence on combining the two integrals in (10), and inverting the order 
of summation and integration, the following expansion is obtained. 


: (1+) 
9 9406 << (° 9 2 
S) (a, B) — (x—f) + mo > A, | 7 : ) -2apsnt df, 
n* " on a ; J (1—?*)! 
Integrating by parts 
; 24 (1+) 
Z 7 "4 r 3 ] ° 
To a S24 _% ______*__|¢-2nysnt d¢, (11) 
7 m2 s J [ie (ep 
Now (Watson, 3, 172) if Rlv } and |argz| < 42, 
K, 7 (32)” * at ¢? -1) + dt 
I(v+4) 
1 
7 2iT'(v+-4) ; 
Hence | e-4(1—#)-* dt = — =~ cos 7K, (2). 


mi(iz 


Analytic continuation extends this formula to all values of v. Putting 


l, 2, in turn, 
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Substituting these results in (11), and substituting the above values of 4 


41,, 


x“ 


9 a 
Ala, 8) = pric —B)- ys (sin 278a—sin 27s) x 


s=1 


] ete tbs 
x | Se128hKy(2mpsn) + P K,(2musn)]. (12 
n 


Substitution of this value for S,(a,8) in (5) establishes the final trans. 
formation given by equations (2), (3), (4) of the introduction. Since 


K (2) ~ (7/2z)te*, as z>+cac 


> 


the series (3), (4), (12) each converge very rapidly provided p is not small. 


In particular : 
Pi , k, ~~ 4mpiste-"7#5, as sax 


These series can be readily computed provided tables of the Bessel 
functions are available. 


3. Numerical tables 

Cases of practical importance include the square tunnel (4 = 1) and 
the duplex tunnel (u = } or 2). For these cases numerical results have 
been obtained; the rélevant values of k, are given in Table 1. The 
computation entailed the subsidiary tabulation of the functions K,(nz) 
and K,(n7). These are given in Table 2, together with the associated 
functions J,(n7) and J,(nz), computed incidentally for checking purposes 
and included for the sake of completeness. With the aid of Table 2, values 
of k, can be computed for all rectangular tunnels for which is an integer 
or half an odd integer. All numbers in Tables 1 and 2 are positive. End 
figures have been correctly rounded off. 

The work described here has been carried out as part of the research pro- 
gramme of the National Physical Laboratory, and this paper is published 
by permission of the Director of the Laboratory. 


TABLE | 





Rs 
s p=4 po 2 p= 2 
I 0°246709 0°026951 0°0001 33 
2 0°013476 0000066 | 0*000000 
3 0:000682 | 0-000000 
4 0°000033 
5 0*000002 


6 ©*000000 





a /6 1'644934 
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TABLE 2 








K,(nz) K,(nz) I,(nz) 1,(n7) 
I 3503 5°47784 50 4°49145 67 
1 0554 )5 6996 87°1085 79°8390 
0003 4216 1°6331 X 10° 1°5439 X 108 
1221 1269 3°261 X 104 3°128X 108 
{ 19 6°7X 10° 6°5X 10° 
© 0002 1°4X 10° 1°4X 107 
O 3X 108 3X 108 
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NOTE ON THE PRECEDING PAPER 
By G. E. H. REUTER 


(Department of Mathematics, The University, Manchester) 


[Received 6 September 1949] 


THE identity found by Olver (1) can be proved by a different method, 
which seems worth pointing out since it could be used in similar problems. 


Olver’s identity is equivalent to 


21—4un > sin 27as{27npsky(2anps)+K,(27nps)], (1) 


1 
vhere f(y) y(y?+- 2u2n?)(y?+-p2n?)-? and > is to be interpreted as 
M gia 
lim > We may write f(y) = n~*f(y/pn), where 
M+0m=—M 
f(a) a(a--+1) Lo (ar*- 1) : 
fx(a~?2+-1) sgn v}+a(#?+-1)-?+sgna 
q(x)-=-sgn 2, say.T 


If we also note that 
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we see that (1) can be rewritten as 
co 
m—x 
g ———— 
S pn 


with 


x 
m 


s=1 


g(x) = {x(x?-1)-++- 





—4un > sin 270s |2amn K,(27nps) + K, (272s) — 


sgn a}+-a(x?+-1)-3, 


To prove the identity (2), we use Poisson’s formula (2) 


Ss h(m) 
where H (x) ( h(t)e 


If we put h(x) = (—), we find that 
pn 


oe 


H (x) = [ ( ‘ —27rixt dt 
pn 


ry 


x 


—C 


2rrizt dt. 


g(v)e—27x( X+ pen “un dv 


— —2 


1.2) 
= —2ipne—27tax { a¢e)sin(2amprn dv, 
0 


since g(v) is an odd function. Now we use (4) and (5) to obtain 


m—o - — are 
g|— = — Zon > e-27ta8(}(2r7uns) 
pn s8=-—a 
m=-—-® 
@ 


= —4un > 


where 
0 


G(x) = | g(v)sin: 
) 


sin 27asG(27pns), 


1 


rv dv. 


To prove (2), it merely remains to show that 


G (2) 


2 @ 


: v ; COS @V 
Now | —sinzvdv = 2x [ <s uae 
} (+i) } (+1) 
0 


so that (6) reduces to 


g 


= he 1}sin av dv = K,(z)—=- 


t See Ref. (3). 


eK y(n) + Kyu) —-, for 2 > 0. 


dv, integrating by parts, 


1 





But 
| ffs 
8 J ((e? 
(2) }° 
(3) 

(4) 


(5) | and sc 





on 
Q tf = 
2. Aa oe 


But 
— 1 \sin xv dv dv sin xv 


| l a . 
a 
nus |’ | \(y2- 1) } 


cos at 


0 





and so (7) is proved 


F. 
s (Cambridge, 1944), 


W. J. OLVER, see preceding paper. 
E. C. TrrcomMarsnH, Fourier Integrals (Oxford, 1937), 


i. 
2. . 
Bessel Functic 


3. G. N. WATSON, 


l 
x | (¢?-+ Ii 


x 


[ dt 
ey 





™ Fs 


> 


sin av dv 


See Ref. (3). 
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AN APPLICATION OF THE METHOD OF STEEPEST 
VOPNTINT a r nm] ~ Tr NT NETO ) 
DESCENTS TO THE SOLUTION OF SYSTEMS OF 
AT T Tyo" T a T Ta T ry TC 
NON-LINEAR SIMULTANEOUS EQUATIONS 
By A. D. BOOTH 
(Electronic Computer Project, Birkbeck College, London) 
[Received 2 November 1948] 
SUMMARY 

The solution of non-linear simultaneous equations using the ‘method of steepest 
descents’ is discussed. A comparison is made with the Southwell and Synge pro. 
cedures for the solution of a set of equations, and it is concluded that the proposed 
method requires only about 1/n times the number of ‘steps’ to reach a solution, 
It is emphasized, however, that in view of the rather complicated expressions which 
must be calculated, it is best suited for use in conjunction with a high-speed electronic 
computer. . 

Finally a method is given for discriminating between real and false solutions in 
the case when the problem is one of minimization rather than of obtaining an exact 
solution. 
lr };(Z,,...; Zx) 0 l,..., M) (1) 
is a set of M simultaneous equations, in the V unknowns (Z,), a ‘solution’ 
will be defined to be any set of (Z,) for which 

M = 

S (9 

ct) ie d; >; 4 

j=1 
is a minimum, ¢; being the complex conjugate of ¢;. This definition is 
convenient in that it covers, not only the case M = N to which the 
classical definition of solution applies, but also the cases M < NV, M>N 
where there may be either an infinity of solutions or no solution at all in 
the classical sense. 

Suppose now that Wee 

$; = F+1f;, 
T ” ‘w 
95 Pj = PGS; (3) 
so that ® is everywhere positive. The problem is now to minimize ® with 
respect to the 2N variables (X,, Y,), where 


Z, = X,4+1Y,. (4) 
For convenience of notation ® will be defined by the relations 
2M : 
ri) 2, pF s (5) 
where: tho; = F(a%4y.0+5 ay), 
bas-1 = Syl ay 6) 
X, = Vor-1 
Y, = Ly. (7) 


(Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 4 (1949)] 
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In order to effect the minimization all that is needed is some systematic 
procedure, such that if any point P,(,...,7,) is given in the 2N dimen- 
sional space defined by the (x,), and if ® is not a minimum at this point, 
then the procedure makes possible the passage to another point P, of the 
space for which ® is less. Repeated application of the process will then 
ead to the minimization of ®. 

It is helpful, at this point, to consider a two-dimensional case of the 
ibove analysis, since here a graphical interpretation is possible which 
renders the argument clearer and illustrates various advantages and 


defects of the different techniques. 











Fie. 1. 


Fig. 1 is a contour map of a part of the function 
@ D(x, y). 


It has minima at P,, Q,, Ry, and cols at C, and C,. Suppose first that P) 
is chosen as the initial point; then to reach the nearby minimum P, two 
procedures are available: 

1, Determine the direction from P, in which ®, considered as a function 
of that axial variable (i.e. y), decreases most rapidly. Proceed in this 
direction until ® is a minimum. In the figure this end-point is represented 
by S,. Having arrived at S,, repeat the process to obtain S,,...,S,. The 
process ey idently converges, in the case drawn, to the minimum F,. 

2. Determine the direction of steepest descent of ®, considered as a 
function of « and y. Proceed along this direction until the function reaches 
aminimum. Iterate until the minimum is reached. Again the process 
converges to P,. These two methods are typical of (1) the Southwell 
relaxation technique (1, 2, 3) and (2) the steepest descent technique (4). 
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The relative merits of the two methods in various typical situations 
will now be examined. In the case of a region, such as that about P, in 
Fig. 1, it is seen that both methods converge fairly rapidly but that, on 
the average, (1) will require about twice as many steps as (2). For a long 
valley, as shown about R,,, the steepest descent from Ry would reach R,, 
in about three steps, whereas the relaxation technique would converge very 
slowly, the approximate drawing suggesting the need for about eighteen 
steps. When a col is present, as at C, or C,, the situation is more complex. 
Kither (1) or (2), as described above, converge to a turning-point of 
rather than a minimum; so that starting from a point Q, in the neighbour- 
hood of C, either (1) or (2) would converge to C, rather than to the true 
minima P,, Q,, or R,,. Similarly, if true maxima of ® are present in the 
region, the processes may converge to these. The above intuitive argu- 
ments suggest that the steepest descent technique is the best for the 
purpose in hand, but that some modification is required to guard against 
ascent to maxima or stabilization in cols. In essence the method, to be 
derived analytically later, consists in making a parabolic approximation 
to the curve of ® considered as (1) a function of (a,) or (2) a function of (n), 
the distance along a normal to® = const., the normal to ® — const. being, 
as is well known from vector analysis, the direction of steepest descent. 
In the first case, 


9 
.; 
O = (0) e,®,+—90,, (8) 
a 
‘ . > can n \ 
and in the second, ® = 0(0)-4 En Pn t+ “On ns (9) 


where ®(0) is the value of ® at the initial point, 
®, = (00/éz,),, (10) 


®,., (¢ *D c x) c Xs)o, II 


€, is an increment in 2,, and e, denotes a change of the coordinates in the 
direction of the normal. For the relaxation version ® has to tend to a 


inimum, whence 
minimum, whence dM /de, 0, 


that is € ®,/®,.,., (12) 


r 


and similarly for the steepest descent 


€ —QD /P (13) 


n n n,n* 

It is easy to see how ascent to a maximum can occur. Suppose that the 
initial point is chosen between the point of inflexion P, and the maximum 6 
(Fig. 2). In this case ®,, is negative as also is ®,,,, (this is equally true for 
the relaxation version); thus (13) gives a negative value for e,, which means 
that P,, moves up the curve towards G. In the reverse situation, when P 
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is at P, between P; and the minimum JL, ®, is again negative, but this 
time ®,,, is positive making ¢, positive so that ® descends towards L. 
\nother fact which becomes obvious from the geometry of Fig. 2 is that 
if P is unfortunately chosen at P,, « will be infinite. These considerations 
suggest that, however excellent the relaxation or descent to a minimum 


techniques may be for homogeneous quadratic forms, where it is known 


C; 















0 M: M, Me MjM)M, (ror n) 


that only minima exist, they are not satisfactory in the initial stages of a 
determination of solutions for non-linear equations of the type given in (1). 

Fortunately another method is available. Approximating ® by a linear 
function of «,,, (9) becomes 


@ (0) MD « ; (14) 


® can be decreased by taking 
‘ @(0)/O,. (15) 


This corresponds to moving down the tangent to the (®,n) curve until it 
intersects the n-axis. 

It is seen from Fig. 2 that for P,, P,, and P, this technique leads to 
ibscissae M,,, M,, and M, at which ® is less than at the initial points. 
It is also seen that if the initial point is too near to the minimum, as at Py, 
the resulting abscissa J, may correspond to an increase in ®. This, 
however, is no disadvantage since a simple working rule is available: 
Use (15) until a refinement is reached which increases D; then change to the 
more accurate formula (13). One other point is worthy of mention in 
connexion with the tangent descent. If the minimum value of ® is actually 
zero as at M,, Fig. 2, and the initial point is sufficiently close to M, for the 
curve to be considered parabolic, then it is easily shown that 


UM, = $M, Mp, 
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so that a more profitable estimate of €,, than (15) is 


€ 20(0)/®,,, (16) 


nt 
and this should be used in preference to the latter. The amount 0(0)—o 
by which ® is changed by the descent processes (12) and (13) is easily 
found by substitution into (8) and (9). For the relaxation method 


0(0)—® tp? 20, »° (17) 


Synge (3) has suggested that instead of taking max(®,) as the direction 
of descent in the relaxation method, a more rapidly convergent process 
would be to take that direction for which the actual descent, given 
by (17), is a maximum. This version of the method involves the deter- 
mination of the ®,,. which may well be a laborious process in complicated 
systems. 

The disadvantage of the full steepest descent method, typified by (13), 
is that, as will be seen in the next section, it necessitates the calculation 
of all derivatives of the type ®,, and requires a considerable amount of 
multiplication. On the other hand, axial descent methods of the type 
suggested by Southwell and Synge tend to become inefficient with a large 
number of variables. It was seen that for two dimensions, about two 
relaxations are, in most favourable cases, equivalent to one steepest 
descent; in N dimensions this number will be of the order N. With an 
electronic computer it would be preferable to use either the tangent 
descent method or the steepest descent to a minimum, since neither of 
these requires much discrimination on the part of the machine. For 
manual computation, however, it would seem that the Southwell relaxa- 
tion process is the least laborious as it involves only the calculation 
of all ®, and one ©... and the discrimination of the largest of the ® 
is easy. 

A simpler method than any of the above, however, is the following: 

1. Determine the values of ® at the point given by applying the 
correction of equation (15), 4(1) say, and at the point given by applying 
half these corrections, 4(4). 

2. Using ordinary quadratic interpolation (or extrapolation if appro- 
priate) on the values 4(0), 4(4), 4(1), find the distance along the normal 
to the minimum. 

3. Repeat process using this point as origin. 

It will be seen later (equation (25)) that the above process requires the 
calculation of only one set of first derivatives and values of ® at three 
points per refinement. 


The normal derivatives required for the steepest descent process will 
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now be evaluated. From generalized vector analysis 


®, = ab/én = |grad| (S @,}', (18) 


COs(”, %,) ®,/®,, (19) 


where (n,x,) is the angle between the normal, n, and the 2,-axis. 


re 


N 


Again ins o(®,,) on = - ,) =" 
Ox, on 
r=1 
2N 2N 
| > %,.0,0,]/[> @,?| (20) 
r,s=1 r=1 
From (18), (20), and (13) it follows that 
2N 3 2N 
i [> @,)]'/[ > ©,,.0,0,|, (21) 
r=1 ! lrs=1 si , 


and consequently, the changes, e,, in individual coordinates (x,), for steepest 


lescent are given by 
. €,, COS(N, X,) 


2N 


[> (®,) Fo,/[3 a , ,®, ®,]. (22) 
r=] 
Similarly the steepest descent along the tangent to the (®,n) curve gives 
/ XY } 
€, o0)/[> (®,)? ic (23) 
from which the changes in individual coordinates are seen to be 
2N 
€. 0(0)0,/[ S (®,)?}. (24) 
I lp=y 


In the interpolative version of the method the values of €, are taken to be 
[(D(1)—4(4)+-30(0)] (0)9,, 
4 D(1)—2(4)+-0(0)| [> (, | 


where (1) is the value of ® at the point given by applying the corrections 


; (25) 


from equation (24), and ®(4) the value with half these corrections. The 


®, and ®, , are obtained from (5), namely 


®. 2> ob; Pb; (26) 

and ®.. 2 > W Bb; Bip thy byes) (27) 

where p;, = Op,;/ex, ) (28) 
Ping = Op, /Ox, Ox, - 


5092.8 


Hh 
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The forms of the %; derivatives depend, of course, on those of the 
particular ¢; of the problem. 

The method has been used by the author in the solution of sets of 
equations of the type 

N. : rd . 38 
4, = [S soonin]'+[5 singe], — 0, 
r=) r=1 
and for three such equations in three unknowns it was found that a solution, 
accurate to 0-01 radian, could be obtained in four steepest descents which 
took about 90 minutes. 

An interesting comparison of the relative efficiency of the various 
methods is to be had via the solution of a set of simple simultaneous 
equations. The equations ry 7 
2aty = 5 
were taken, and assuming initially that « = y = 0, it was found that 
eighteen applications of the Southwell or Synge relaxation technique gave 
the values 2 = 0-98. y = 3-02, ® — 0-0008. 

Five applications of the tangent descent method gave 
2 = 1-04, y = 2-97, ® = 0-003. 

After which ® increased slightly, indicating breakdown of the method in 
the manner of My of Fig. 2, and finally, four applications of descent to 
a minimum or the interpolative method gave the correct solution: 

x 1-00, y 3-00, ® — 0. 
During the calculation of the relaxation version it became obvious that 

> 

the convergence was very slow in the region of the minimum. 

The five procedures for obtaining a solution of (1) via minimization 
of (2) may now be summarized. 

1. The Southwell relaxation method. Find that axis (x,) in the 2N- 
dimensional space of the variables (x,) for which ®, is greatest. Change 


x, to (x,+e,) where 


€ —,/®,., 


; 
and repeat the process until ® is minimized. 

2. The Southwell—Synge relaxation method. Find that axis (x 
07/20,,. is greatest. Change x, to (x,+e,) where « 
Repeat until ® is minimized. 


) for which 


a 
, is as given in l. 
3. The steepest descent along a tangent. Change all coordinates (x,) to 
(x,+-e,) where the e, are given by 
2N 


, = —0(0)0,/[ > @,?]. 


Repeat as in 1 and 2. 
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1. The steepest descent toa minimum. As in 3, but with the e, given by 


é, [> (©, )*|®, / | S ; ©, ,0,0,]. 


r,s 


5. The interpolative descent. Calculate corrections as in 3. Evaluate ® 


at the point given by applying these corrections (0(1), say) and also at 
the point given by applying one-half of these (0(4)). The corrections are 


then taken as £ [(1)—4(4) 4 0)] (00, 


4/O(1) 20(4) 4 “hi ))| [> 3] 
One point remains to be mentioned, the desirability of knowing whether 
the method has converged to a col, as at C, and C, in Fig. 1. In the case 
of a system having solutions in the classical sense, convergence to a col, 
rto a relative minimum, can be seen directly from the fact that ® is not 
zero. When, however, the problem is one of minimization only, the 
distinction between cols and minima becomes important. The simplest 
test is to find ®(min) and then to consider the quadratic form 
2N 
a) > €, e,®,, ,. (29) 
r,s=1 
If ®(min) is a true minimum, then Q is a positive definite form. This fact 
can be ascertained by the following procedure. In matrix notation 
OG = «Ae, (30) 
where e’ is the transpose of « and A is the matrix of the _ 
Let C be an orthogonal matrix, then a transformation 
«= Ce (31) 
can be found such that 
Q = «'Ace = e'C'ACe = e'Ke, (32) 
where A is a diagonal matrix, the quadratic form being thus reduced to 
2N 
Q= > L,,é. (33) 


a 


r=1 
Now let k be any element of K; it is evident that 


K—kI| = 0, (34) 
whence, irom (32), C’AC—kI1| = 0. 
But since C 1s orthogonal c'¢ CU : l, hence 


C'’AC—kI| = |C'AC—C'kIC| 
C'(A—kI)C| 
= |C’|.|A—kI|.|C| 
- |A—Kkl1| 
A- = 0. 


whence 


kI| 
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This shows that the elements, k,,, of the diagonal matrix K are the latent 

roots of A, i.e. the solutions of the secular equation 
| @y,—k Ay Q1 on 


| 
| 
| 
| 


a,,.—k ‘ = 0, (36) 








Mon1 Gav on—k 


and, from (33), it is evident that Q is positive definite if and only if all 

kv : 0. 
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THE CONDITION OF CERTAIN MATRICES. 


By JOHN TODD 
(National Bureau of Standards, Washington, 25, D.C.) 
[Received 3 May 1949] 


SUMMARY 
Two measures of the condition of a certain matrix are calculated. The matrix is 
that of the system of linear algebraic equations arising in the approximate solution 
of the partial differential equation 


| | i 
[as 2s kz. 
ox” oy” 
The results obtained confirm analytically the experimental fact (observed by com- 
puters) that the condition of this matrix is better than that of the corresponding 
matrix in the case of the ordinary differential equation y” = ky. 


1. THE condition of the matrices associated with the numerical solution 
of ordinary differential equations of the second and of the fourth order has 
been discussed recently (1), (2). We now consider the condition of the 
matrix associated in a similar way with the solution of a second-order 
partial differential equation of elliptic type. For simplicity we consider 
the solution of 

—to5 = (1.1) 
ox" oy” 

within a square, of side /, on the boundary of which z vanishes. 

It will appear that the solution of the simultaneous linear algebraic 
equations involved in the approximate solution of this equation will be 
relatively less troublesome than the solution of the corresponding system 
in the case of an ordinary differential equation. In fact the first system 
is not ill-conditioned while the second is slightly ill-conditioned (1). 

2. We divide the square, the vertices of which we assume to be (0, 0), 
(0,1), (7,0), and (7,2), into (n-+-1)* squares each of side h, by the lines x = rh, 
y = sh, where r = 1, 2,..., m and s = 1, 2,..., n. If we replace the partial 
derivatives by a finite difference approximation, e.g. @z/@x? by 

h-*{z(x-+-h, y) —22(x, y)+2(4—h, y)}, 
we obtain instead of (1.1) a system of n? equations in the n? variables 
2,3 = 2(rh, sh). The general equation of this system is 
Zpmt,s 1 2p,s—1— pg 1 2418+ 2841 = h*rz,,g 
(there are slight modifications at the boundaries, i.e. when r or s is 1 or n). 


If we now put z,, = 2,.»(,-3) and denote by Z the column vector whose 


[Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 4 (1949)] 
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j++, 2, the system of equations can 


— are ” 7 9 
1 n Somponents are the Z4,~=1,2 
be exhibited in the form 











AZ=M’*Z or, say, BZ=0, 
a Sr ae a: (_—4 1 
| a a Geer eee ) =~ ] 
a oa ae : 1 4: J ; 
o.oo ; X : : i 
gers 2 2 ee ~« « Sabo 
a see eros ae aoe ; i © « l = 


Here X is 


matrices. 


an nxn matrix and A is a vv matrix, partitioned into nx» 
Tisann Xn unit matrix. 
We shall discuss first the condition of the matrix A and later show that 


the condition of B is not essentially different. 


3. It follows from investigations of von Neumann and Goldstine (3) 
that a reasonable measure of the condition of a matrix is the ratio |A\/\y 
where A and wp are respectively the greatest and the least (in absolute 
value) of the characteristic roots. We have called this ratio the P-condition 
number of the matrix.: When it is large the matrix is ill-conditioned. 

It has been shown (4) that the characteristic roots of A are 


r..g = —4—2 cos r6—2 cos 80, 
where @ = z/(n+1) and r = 1, 2...., n, 8 = 1, 2,...,n. It follows that the 
P-condition number of A is 
4+-4cos 6 8 4n2 


O(v). 


> y — lt > A —_ >» 
4—4cos@ 8 sin? 16 7 
Thus the condition number of A is of the same order as that of an average 
matrix (see (3), p. 1097). Further, A 
corresponding matrix in the case of an ordinary differential equation 


is in better condition than the 


which has a P-condition number of order v2. 


4. Turing (5) has suggested the N-condition number, defined as 


N(A)N(A-)/, 


oe ro]? . ee : 
where N(A) = [> > as|°| , as a measure of the condition of a vx» 
1 


r=1s= 
matrix A. This number can be estimated quite accurately in the present 
case. It is clear that 


N(A) = (20n?—4n)? ~ V20n. 








The explicit inversion of A and the subsequent direct computation of 
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N(A-1) seems best avoided and we use an indirect method. We recall 
» fact that — ee 

the N (¢ ) [> A? }t. 

where C is any matrix and the summation is over all its characteristic 

roots A (see e.g. (3), p. 1046). We also use the fact that the characteristic 

roots of A-! are the reciprocals of those of A (see e.g. (6), p. 23). It 


follows that, in the notation of § 3, we have 


N(A-") [> >a] 
r=1s=1 


. ne nr ae = ol 
whence N(A-) | > > {sin’r¢ -sin?sp}*] : 
r 8 1 
where d = 30 = /2(n+1). 


It is well known that 
x >sine > 2a/r 
provided that 0 < # < 3m. Hence, each summation running from | to n, 
[> > (r?+8?)-24-4]# < 4N(A-1) < Jr ¥ Y (r?-+82) 2b]. 


The double sums can be estimated in a familiar way: first replace them 
by double integrals with respect to x, y over a square and then by double 
integrals with respect to r, 6 over quadrants included in and including 
the square. We find that N(A-) lies between two fixed multiples of n. 
This means that the NV-condition number is O(1) exactly, which is smaller 
than that of a random matrix which is O(v+) (see (5), p. 299). Again, 
there is an improvement on the corresponding results for an ordinary 


differential equation which was O(v?). 


5. It remains to discuss the actual matrix B. We observe that 


We assume that A is a characteristic value of the system. The character- 
istie roots of B are obtained by subtracting Ah? from those of A which 
ire the A, of paragraph 3. This subtraction will not influence the order 
of the P-condition number but only alter the constant implied in it. We 
have, in fact, for h sufficiently small, that the P-condition number of B is 


4-1-4 cos @—Al?n-* 8n* 
m “ od aes : ‘ : O(v). 
8 sin?é—Al?n-* 27? — Al? 
Neither will the order of the N-condition number change for we still 
will have V(A) ~ V¥20n while the subtraction of Ak-? will involve changes 
in V(A-!) only to the extent of changes in the constants between which 


n 1N(A 1) lies. 
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6. The cases when an equilateral triangular net or a hexagonal net js 
used instead of the square net can be investigated in a similar manner, TH 
Considerations similar to the present can be applied to the more general 
case m. 22, 
r=1 " 


It is found that as m, the number of variables, increases, the condition of 


i 
the corresponding matrix improves; as m — oo the condition approaches 
that of the unit matrix, which is perfect. Detailed discussions of these 
matters will appear elsewhere. 

Fi 
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THE NUMERICAL METHOD OF CHARACTERISTICS 
FOR SUPERSONIC FLOWS WITH AXIAL 
SYMMETRY 
By M. HOLT 


(Department of Applied Mathematics, The University, Liverpool) 
[Received 8 March 1949] 


SUMMARY 
If excessive labour is avoided, the numerical method of characteristics, in its 
present form, gives less accurate computed results in fields of flow with axial sym- 
metry than in fields of plane flow. In this paper a simple process of error correction 
is developed by which the accuracy in fields of the former type may be improved. 


1. Introduction 

We refer a field of steady, supersonic, isentropic, irrotational flow, involv- 
ing two space coordinates, to rectangular Cartesian axes (with the x-axis 
coincident with the axis of symmetry for axially symmetrical flow) and 
choose the Mach angle, », and the angle between the direction of the 
velocity and the x-axis, 6, as independent variables. The governing 
equations of the flow, expressed in characteristic form, are then (see, for 


example, ref. 1) dy dx = tan(@+ 1), (1) 
f(u) du db + jg*(9, n) dyly = 0, (2) 


where f() = 2cos*u/(y—1+-2sin’x), g*+(0,) = sin @sin p/sin(@Fp), and 
j = 0 in plane flow, 7 = 1 in flow with axial symmetry. 

The numerical method of characteristics is one by which equations (1) 
and (2) may be integrated step by step. A given non-characteristic initial 
line is firstly divided into a number of small segments. At the vertex of 
the characteristic triangle based on each segment approximate values of 
the coordinates and of the physical variables are found from the initial 
data, using equations (1) and (2). In this way a second line is constructed, 
from which the process may be continued until the state of motion in the 
whole region determined by the given line has been found. 

Figure 1 shows an element of the characteristic network in this region 
of influence. Points 1 and 2 bound any segment of an initial line, and A 
is the point of intersection of the plus characteristic through 1 and of the 
minus characteristic through 2. The application of the method to plane 
fields on the one hand, and to fields with axial symmetry on the other 
hand, differs in the process by which the state of motion at A is determined 
approximately from the known state at points 1 and 2. 


[Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 4 (1949)] 
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In plane flow, owing to the absence of its third term, equation (2) yields 
universal relations between 6 and yp along characteristics. In consequence, 
if the values of 6 and » are known at one point on each of two charac- 
teristics of different families, the values at their point of intersection can 
be found exactly. The difference process employed in plane flow, due to 
Busemann and Prandtl, and set out by Temple (ref. 2), takes advantage 


1 





2 A 


Fie. 1. 


of this property. With the values of @ and » at A known exactly in 
advance, the coordinates of A are found approximately by taking the 
slopes of the chords 14 and 2A to be the means of slopes at 1 and A and 
at 2 and A, respectively. 

In flow with axial symmetry it is not possible to separate the determina- 
tion of 6 and pw from that of # and y; in fact before the characteristic 
directions at A can be found, the y-coordinate at A must be known. The 
use of a mean difference process is thus excluded, at least in the first 
approximation. Instead the state of motion at A is found by a simple 
difference process; the coordinates of A are found as those of the point 
of intersection of tangents to characteristics at points 1 and 2; the y- 
difference along each characteristic is then found, and the values of 6 and u 
are determined from equation (2), replaced by a simple difference equation. 

The second process described is less accurate than the first. It may be 
shown (cf. the appendix) that, if the length of the initial segment 12 is 
of order e, the first process introduces errors in 2, y, 6, and uw at A of order 
e°, and the second process introduces errors of order e?. 

In principle, the accuracy of the second process could be raised by 
iteration, that is, by performing a second integration using mean differ- 
ences calculated from the results of the first integration. The results set 
out below give support to a different approach in which errors of order ¢ 
may be either calculated or eliminated with the use of a simple difference 
process alone. 
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2. Errors in x, y, 0, and » at the vertex of an elementary charac- 
teristic triangle 
We now give formulae, correct to O(e?), for the coordinates and flow 
variables at the vertex of an elementary characteristic triangle in an axially 
symmetrical field. They involve the difference between two successive 


1 


A 
Fic. 2. 
approximations to conditions at the vertex, the first starting from a single 
interval on the initial line, as described in § 1, and the second starting 
from two half-segments of the initial line. 

In Fig. 2 the points A, 5, and 6 are the vertices of triangles based on 
initial segments 21, 01, and 20 respectively. In the first approximation A 
is reached by a one-step process from the single initial segment 21. The 
second approximation employs a two-step process; starting from the two 
half-segments 20 and 01 of the initial line conditions at the corresponding 
vertices 6 and 5 are found by the one-step process, then the next approxi- 
mation to A (not shown in the figure) is reached in one step from 6 and 5. 

The position in the field at which a variable is evaluated is defined by 
an appropriate suffix. The order of the step process used to estimate the 
variable is indicated by an index. Values found correct to O(e*) are 
written without an index. The slope of a plus characteristic, tan(@—,), 
is denoted by S; that of a minus characteristic is denoted by 7’. Finally 
0 is chosen so that x,»—2, = %,—%_ = he. 

With this notation, it is shown in the appendix that the correct values 
of x, y, w, and @ at A are 


ty = e4+2(x2,—2z}), (3) 
Ys = Yat2(y4—Y's), (4) 
M4 py +2(u2,—py), (5) 


0, = 0+ 2(62,—64). (6) 
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It is further shown that if 


U = i(x4—2,)(S4—S)), (7) 
V = }(x4,—2,)(T4—T,), (8) 
L = —}(f4—f,) Meat) 1HGa /Ya—9t lYYa—-Y), (9) 
—4(f4—f2)(M4— M2) 49920 |Ya—92 [Yo)(Yu—Yo), (10) 

then 
t4 = t4+2(U—V)/(T,—S,), (11) 
Y4 = Ya4t+2(UT,—VS,)/(T,—S,), (12) 


pet = (L+M)/(fi-fe)—ILGF ys) + Ge /y2) ¥4—y4)/(fitfe)], (13) 
P—F = (M—L)/(f,+f2)—I(fe9t /Ys)—(fi 92 /¥2) I (ys—ya)/ (fr +F)].- 
(14) 

These results show that we may calculate errors in x, y, p, and @ in 
either of two ways. We can find the second approximation to values of 
x, y, w, and @ at the vertex by a two-step process and then employ 
equations (3), (4), (5), and (6). Alternatively, we can use equations (11), 
(12), (13), and (14) to express the errors solely in terms of the initial data 
and the data obtained by the one-step process, and so avoid an extension 
of the difference process. 

If the errors at the vertex of a triangle based on one initial segment are 
known we can easily find errors at vertices corresponding to adjacent 
segments of the initial line. In fact, if A and B are the vertices of two 
adjacent characteristic triangles based on segments of the initial line of 
length ¢ and 6 respectively, it can be shown that, for each variable 2, y, 


H, and 6 (error at A)/(error at B) = e2/82. (15) 


It follows that the full process of error calculation need only be applied 
to a limited number of the triangles based on the initial line, errors in the 
remaining triangles being found by use of the simple relation (15). The 
actual number of triangles requiring use of the full process will depend on 
conditions on the initial line, in particular on the derivatives of y and é 
along the line. 

The author wishes to thank Professor 8. Goldstein for his interest in 
this paper; this forms part of the author’s Thesis for the Ph.D. Degree 
of the University of Manchester. 
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SUPERSONIC FLOWS WITH AXIAL SYMMETRY 
APPENDIX 

To establish the error laws stated in § 2 we need to generalize the process of 

approximation and to find conditions at the vertex A of the triangle in Fig. 2. 

The second approximation begins with the calculation of conditions at 6 and 5 in 

one step from 20 and 01 respectively; in the third approximation conditions at 6 

and 5 are firstly found in two steps from 20 and 01. In general the nth approxima- 


tion starts from 2”~! intervals of the initial segment; conditions at 6 and 5 are found 


in 2”-2 steps from 20 and 01 respectively, and we then proceed to A in 2"—* steps 
from 65. 

Firstly, we prove relations (3) and (4). The correct values of the coordinates of A 
are the limits of 2 and y", as n — oo. To relate these limits with data obtained by 
the more simple processes, differences along plus characteristics and minus charac- 
teristics are considered separately; firstly we treat the line 15A. 


First approximation 
From equation (1) and the initial data we find 


yy —Y, (a, —2,)S,. (A. 1) 


Second approa imation 
We use the first approximation to give differences between 1 and 5 and then 
between 5 and A 


Ys—Y1 = (x}—2,)S), (A. 2) 

Yu —Ys = (74 —275)S}. (A. 3) 

Hence Y4—Yy (x?, —x,)S,+ (a3, —2})(S$—S)). (A. 4) 
In equation (A. 4) only terms of order e? are retained. Each factor of the second 


term on the right is of order e; in evaluating these we need therefore retain only 


terms of order e, and employ only the first approximation. Then, correct to O(e), 


e2, — 2}, al, — ai 4 (a —z,); 
S; S; 3(S%—S)), 
and, with U as in (7), equation (A. 4) becomes 
Y4— Yr = (%4—%,)S,+U. (A. 5) 
Subtract (A.1) from (A.5). Thus 
4A Y's (4- x})S,+U. (A. 6) 


In equation (A. 6) y2,—y, and x*,—2}, are both of order e*, so the factor S may be 
evaluated at any point in the characteristic triangle. To be consistent with later 
results we shall use its value at 1. 


If we introduce ’ 7 
e introduc D yk —ak S,, (A. 7) 


= U. (A. 8) 


equation (A. 6) becomes D*, _ D} 


The nth appre rimation 


It can be proved by induction (using arguments similar to those given above) that 


the equation corresponding to (A. 6) in the kth approximation is 
j » 6 rs k-27 17 
Ya— Yi (ark x,)S,+[2—($) ]U. (A. 9) 
Putting } n, n—1 successively in (A.9), and taking the difference between the 


resulting equations, we find that 


Di, — Di-* = (4)"*U. (A. 10) 
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Adding equations (A 10) when n takes all positive integral values from 2 to n, w, 


obtain : ; m : a 

Di — D3, = (1+(4) + (92+. + (410. 
Hence, D, = D,+2U. (A. 1] 
Similarly, if E" yi—x, T, we can show, by considering differences along thy 


minus characteristic 26A, that 
E, = E4+2V. (A. 12 
Equations (11) and (12) follow from (A. 11) and (A.12). Equations (3) and (4) an 
obtained by combining (A. 11) and (A. 12) with (A. 8) and the corresponding equation, 


EH? — EH}, = V. (A. 13 


Equations (5), (6), (13), and (14) are established by repeating the above process of 


approximation starting from equation (2) instead of (1). 
Note 

We have assumed throughout that the error introduced during each step in th 

g 
process for axially symmetrical flow is of order e?. This result may be proved by 
repeating the analysis used to establish equations (A. 8) and (A. 13) retaining only 
terms of order e; we then find 
Di{—D,=0, E,{—E,=0, sothat z,=2% and yl = y2. 

In similar manner, we can show that the error introduced at each step by the mean 
difference process for plane flow is of order e*. 
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THE PENDULUM 
By L. M. MILNE-THOMSON (Royal Naval College, Greenwich) 
Received 31 August 1949] 


SUMMARY 


The discussion of the simple pendulum by means of elliptic functions is an old 
problem. Nevertheless the following concise treatment may be of interest. 


THE energy equation of a pendulum, simple or compound, can be expressed 


in the form 
G2 w*(1—k* sin? 36), 4q/(lw?), (1) 


where @ is the inclination to the downward vertical, w is the value of 6 
when @ = 0, and / is the length of the equivalent simple pendulum. Take 
i as the modulus of Jacobian elliptic functions and write 
sin 3}@ = snu = sn(u,k). (2) 
Then cos 36 cn u, 6 = wdnu. (3) 
Differentiate (2); since d(sn u)/dt = cnudnuv we get, with the aid of (3), 
i = ha, or u Lwt if 9 = 0 when ¢t = 0. Therefore 
sin 46 sn{ Sut, k), k? 4q (lw?) (4) 
which solves the problem of the pendulum in terms of elliptic functions. 
Case (i) The pendulum makes a complete revolution in time 7’. In this 
case 6 increases from 0 to 7 while ¢ increases from 0 to $7’ and therefore 
sn{ wT’) a bw T' = K, 9 4K Ww, 
where K is the real quarter period of the elliptic functions. 


Case (ii) The pendulum oscillates through the angle « on each side of 


the vertical. Here 0 0 when @ x and so 

2 > 102 1 1 »Qsec 1 / 

/ cosec* 5a, sw COSEC 5a J (g/L), 
while from (4) sin 36 sn(twt, cosec 4a). 


Since the modulus now exceeds unity, we use Jacobi’s real transformation 


to the reciprocal modulus to give 


sin 30 = sin }asn(}wt cosec $a, sin $a) 


, gq . 
sindasn{t /2,sin da}. 
7 ( Ji , 5 ) 


(Quart. Journ. Mech. and Applied Math., Vol. II, Pt. 4 (1949)] 
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If T is the period of an oscillation, @ increases from 0 to « while ¢ increases 
from 0 to 47’. Therefore 


afin f)—1, roan ft 
g 


where K is the real quarter period corresponding to the modulus sin 4y, 
The value of K is obtained in either case from tables of complete elliptic 
integrals. 
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